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INTRODUCTION 


Recent satellite-borne radar altimeters such as GEOS-3 and SEASAT-1 have Included 
waveform sampling gates providing point samples of the transmitted radar pulse after It 
has been scattered froir the ocean's surface. There were 16 such samplers In GEOS^J and 
63 on SEASAT-1, and Table I lists the time locations of the 63 SEASAT-1 samplers relative 
to the middle one of the set. The waveform sampler data can be averaged over some speci- 
fied number of individual radar returns, and the resulting set of average values can be 
best-fitted in some sense by varying the parameters In a model or theoretical mean return 
waveform. This report describes the waveform model used, and the details of fitting this 
model to the waveform data to obtain estimates of the parameters characterizing the 
modeled waveform. 

The particular model radar return waveform used at Wallops Flight Center for SEASAT 
data analysis is characterized by six parameters: 1) amplitude; 2) time-origin, or track- 

point; 3) ocean surface rms roughness; 4) noise baseline; 5) ocean surface skewness; and 
6) attitude angle, or off-nadir angle. Figure 1 sketches the model waveform and Indicates 
qualitatively how the model shape changes as these parameters vary. These will be dis- 
cussed later, but we can briefly state here the Importance of these estimated parameters. 

The time-origin parameter is the location of the actual mean return waveform relative 
to the altimete'^'s altitude-tracker-positioned sampling gate set, and the time-origin Is 
thus directly interpretable as an altitude correction to be applied to the real-time 
altitude output, 'he ocean surface rms roughness provides a revised estimate of the 
significant waveheight (SWH). The ocean surface skewness parameter provides additional 
information about the surface elevation probability density function and possibly also 
about the ocean wave spectrum. The amplitude parameter may be used to revise the altim- 
eter-estimated surface backscattering cross-section o®, and the attitude angle also leads 
to a correction to The noise baseline parameter Is of relatively little direct 
Interest but must be Included as one of the fitted waveform parameters because the sam- 
pling gates measure radar signal plus noise. 


TABLE I. TIME LOCATION AND INDEXING FOR THE 63 SEASAT WAVEFORM SAMPLERS 
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Mean Return Waveform 


The S£ASAT fitted waveform Is characterized by the following 6 parameters: 
Amplitude. 

Time origin. 

Ocean surface rms elevation o . 

Baseline. 

Oc^an surface skewness, and 
Off-nadir angle 



Figure 1. SEASAT model waveform. 


This report Is Intended to be relatively complete and self-contained In describing 
the waveform model and the parameter-estimating procedure. The work depends of course on 
a number of other reports and papers which will be referred to In the texti but all neces- 
sary mathematics should be contained In the report or In the program source listings In 
the Appendix sections. The actual computer Implementation Is relatively simple but It Is 
not S',.‘ simple to describe succinctly the number of judgments and decisions built Into the 
work. I have tried to refer to some of these considerations In this report and this Is my 
partial excuse for a more wordy account than tiiat of a typical technical journal article. 

There are several separate categories which would be Included In a oxnplete account 
of our SEASAT data analysis. Including the following topics: 

1) details of waveform data Input to the programs, data calibration and correction, 
data editing or screening, the waveform averaging period, and so forth, 

2) the general procedure used In fitting a several -parameter function to the 
corrected waveform sampler averages, 

3) the details of and justification for the model waveform used as the several- 
parameter function In the fitting procedure, 

4) verification of the model and of the fitting programs, and 

5) presentation of specific results and ground truth comparisons for the results of 
the data analysis. 

This report Is tm a complete SEASAT data study and will not present apy of the 
Information In 5) above. Item 1) will be Ignored except to note that most of our SEASAT 
work typically used 10,000 pulse, ten second waveform data averages. Item 4) on verifica- 
tion Is only briefly discussed later. This report will concentrate on Items 2) and 3) In 
the following two sections, describing first the general Iterative non-linear least- 
squares fitting procedure and then the model waveform used. At the heart of any waveform 
data fitting procedure Is the theoretical model which Is being fitted, and this report 
places Its major emphasis on that topic. 

As another way of establishing the Intended purpose of this report. Figure 2 shows 
the various procedures which should be Included within the general category of waveform 
processing In planning the data processing for any future satellite-borne radar altimeter. 
Figure 2 Is a rough sketch based on some current Wallops Flight Center studies for the 
proposed TOPEX radar altimeter program, but the figure Is general and not restricted to 
any one specific program. Outputs from the general waveform processor of Figure 2 In- 
clude; 1) new quantities not otherwise available (such as ocean surface skewness); 2) 
quantities which are replacements for, and are more precise than, radar altimeter real- 
time estimates (such as SUH); and 3) quantities which are corrections to altimeter real- 
time estimates (such as altitudes). Within the general waveform processor of Figure 2 Is 
a box labeled Waveform Modeled Parameter Recovery , and this present report Is concerned 
only with this one particular box In the general Figure 2. 
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Figure 2. Procedures within overall waveform processor. 

There are a couple of points of general philosophy to mention here. First, we have 
strongly preferred fitting the data from all of the SEASAT waveform samplers in the 
plateau region of the return pulse. Otherwise one would need to decide what samplers to 
ignore, and this criterion would probably have to be a function of significant waveheight. 
Second, the model waveform is the convolution of several terms as described later, and the 
parameters of interest could be obtained by some type of deconvolution process as done by 
Lipa and Barrick (ref. 1). A direct (nunerical) deconvolution in this type of problem, a 
problem characterized by an inherently noisy set of discrete observations over a finite 
time span, has the general reputation of significant numerical difficulties and insta- 
bilities. Accordingly, at the onset of this work some years ago, we chose to use what 
could be characterized as "iterative convolution," in which an assumed functional form is 
fitted through a non-linear least-squares procedure, and this is the method described in 
the remainder of this report. 
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Finally, a brief paragraph on the history of the SEASAT analysis and on the relation- 
ship of this report to othe'* work may be useful. There has been an active effort In 
waveform analysis at Wallops Flight Center since before the 6E0S-3 satellite was launched. 
The waveform analysis used here for GEOS-3 was described briefly In an earlier report of 
mine (ref. 2), and this algorithm Is among several which were compared In a paper by Fedor 
e^^. (ref. 3). Our earlier GEOS-3 work and the present SEASAT-1 work are based on the 
radar-ocean surface Interaction sunmarlzed by Brown (ref. 4). In a recent paper (ref. 5), 
I developed a several-term expansion describing the mean return waveform from an Idealized 
SEASAT-llke radar altimeter; the paper's results do not apply directly to the actual 
SEASAT-1 radar altimeter because the actual SEASAT-1 radar pulse (sampled by a calibration 
mode) Is not the simple form assumed In Reference 5. Consequences of the actual SEASAT-1 
point target response and of certain calibration problems are described In a report 
(ref. 6), and preliminary results of our waveform processing of SEASAT-1 data are In 
References 7 and 8. Details of the actual SEASAT-1 altimeter design and altimeter hard- 
ware are available In References 9 and 10, and some of the standard ground-based data 
processing for SEASAT-1 data Is described by Hancock et al,. (ref. 11). A later section of 
this report discusses numerical convolutions performed through use of the Fast Fourier 
Transform, and Brigham (ref. 12) Is one of many sources of information on this subject. 

We have already mentioned the somewhat different approach to waveform processing by Lipa 
and Barrick (ref. 1). The principal goal In our work has been the estimation of sea 
surface skewness from satellite altimeter data, and the work of Huang and colleagues 
(refs. 13 and 14) provides a relationship between ocean surface spectrum details and 
skewness and significant wavehelght. 


NON-LINEAR LEAST-SQUARES FUNCTION-FITTING ROUTINE 


Introductory Discussion 

The problem of fitting a function to Input experimental data has been divided Into 
two pieces; this section of the report will describe the fitting procedure Itself with 
little concern for what function Is being fitted, while a later section of the report will 

describe the specific function (the model waveform) to be fitted In the radar altimeter 

waveform problem. 

A least-squares fitting procedure Is desired but there Is a limited number of ty s 
of functions for which the normal equations are directly solvable In an ordinary least- 
squares procedure. The model waveform function described later In Utls report produces an 

Intractable set of normal equations, and In this situation the approach Is to: a) make a 
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first guess at the values of the parameters In the fitting function* b) do a first-order 
Taylor series expansion In parameter space around the first guess values* c) do a least- 
squares solution for the corrections to the first-guess parameter values* and d) use the 
revised parameter values as a new first guess* looping back to step b) and repeating until 
some exit criterion Is satisfied, litis was used earlier In the GEOS-3 data analysis 
(ref. 2) and the general approach has been described In a number of different sources such 
as Reference 15. (References 16 and 17 provide some additional discussion of some Inter- 
est to this general problem.) The following subsection will present this method In more 
specific detail and a subsequent subsection then discusses the particular subroutine 
version Implemented here. 


General Derivation 

Suppose that there Is a set of experimental data Y^* 1 « 1 to N* with each from 
a different value of an Independent variable t^. Suppose further that a model function R 
Is to be fitted to these data, where R Is described by a set of constants K and a set of m 
parameters 0 to be fitted* 


0 = {e,. 02. . . . e^} . (1) 

We will designate the model function value at t^ by R(t^; K* e). The number of points N 
must be greater than the number of parameters being fitted* m. For SEASAT* R will be the 
model return waveform* the Y^ will be the N « 63 waveform sampler values, and number of 
parameters fitted will be m ° 6, but the derivation In this section Is Intended to be 
general and not specific to SEASAT alone. 

Let 0^ be the set of Initial cuesses for the m parameters, 

©0 ' { 0 , . 02 ... . ) . ( 2 ) 

0 ‘‘0 0 

and expand R In a Taylor series about the point In the m-dimenslonal fit parameter 
space. For notation convenience, define 

R. » R (t.; K* e ) , (3) 

0 

and the Taylor series Is given by 
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9R. 


3R. 


(*i* ®> * "i- * ) 15 “ * <*2-*2,> isr 

0 0 1 0 c 


dR^ 

... + ♦ higher-order terms. 

0 m 


Define j and Sj by 


3R, 


’u ■ 55]" * 55J *<*i* ®> 


evaluated at 0. 


( 4 ) 


(5) 


and 


‘j ■ • <«) 

Ignore all the higher-order terms, and the above Taylor series becomes 

R(t^l K. 0) ■ R^ + ♦ ^2^12 ^ ^ ^m^lni * 

0 

This Is a linear equation In 6^. and can be treated by a least-squares method to 

solve for the 6i...6„. 

Although most of our SEASAT analysis has been for uniformly weighted data we will 
Include the possibility of nonuniform data weighting In the following. Define a set of 
weights w^ for the (Y^. t^) experimental data points, and then define Q as the weighted 
sum of squares of the residuals between the and the R(t^; K. 0 ). 

Q ■ I '"^l ^^1 ~ ®)3 

1«1 

- Z w^ [Yi‘^^-^iD<i“«2Di 2 - • • • • Vlm^^ * 

A minimum In Q 1$ found by taking the partial derivatives of Q with respect to *** 

and setting each of these to zero. This set of m partial derivatives Is 


- ... - 6^D.^) - 0 
1^ •■ - 2Ew^ ^l®n“^2®1? * • • • - ® 

( 9 ) 

• • • 

• • . 

1^ » - 2Ew^ ^l‘^n‘^2^12 - • • • * Vlm^ * ® * 

All the Indicated sums are from 1 « 1 to N. These m equations can then be written as a 
matrix equation. 


B • A A 


( 10 ) 


In which B and A are m by 1 column matrices and A Is a symnetrlc m by ra matrix. These 
matrices are 


SWiOnCY^-Ri ) 


and 


m 




o_ 


(11) 


A ■ 


EWiDn' 


ZWiDi,0i2 


EWi0i,D^2 ^1^12* 


^‘*l‘^11®1m 


^N^ll'^lm ^’^l^^^lm * * ‘ ^“l® 


1m 


( 12 ) 
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03) 


Then If matrix A has an Inverse A“\ one solves for A‘^ and the solution for A Is 

A ■ A"^ B. 

Then from the derinitlon of the 6j In equation (€)• the new parameter estimates 6j can be 

obtained from 6. and by 
j J© 

6j ■ 0 j ■»“ 6j , J ■ 1 to m . (14) 

These new 6j parameters are not necessarily equal to the "true” or final paruneter set 0 
because of the higher order terms having been Ignored In the Taylor series leading to 
equation (7). The next step Is to use the of equation (14) as a revised guess» a 
revised 9^, and to repeat the entire process Just described. 

This Iterative process Is repeated for some preset maximum number of Iterations or 
until some type of convergence criterion Is satisfied. Additional error exit provisions 
must be made for cases In which the matrix A Is nearly singular or in which Q, the 
weighted sum of residuals squared. Increases over Its value In the preceding Iteration. 

The usual convergence criterion Is that the fractional change In Q from one Iteration to 
the next should be less than some preset limit. An additional Iteration exit should be 
provided for the absolute value of Q falling below some lower limit even If the Q frac* 
tional change criterion has not been met. The particular lower limit specified Is based 
on reasonable Judgment In the particular application of the general fitting procedure, and 
also can be Influenced by computer precision limits (hje to finite computer word lengths. 


Additional Discussion of the Wallops Version of 
a Non*1 Inear Least- Squares Subroutine 

The subroutine SWHFIT In Appendix B carries out the general Iterative least-squares 
procedure described above. These several paragraphs will d1sv.uss some of the features and 
considerations In this subroutl.ie. 

The function being fitted is provided by a subroutine FILLV udilch Is described 
further in this report's section on the model waveform. The various partial derivatives 
with respect to the fitted parameters, the D^j of equation (5), are supplied by a sub- 
routine FILLO. Notice that by keeping FILLV and FILLD outside the matrix-fitting loop 
within SWHFIT, we have the ability to modify the model waveform function at will without 
disturbing the fitting procedure; also the FIUV function Is evaluated for all N (usually 
63) values at one point In SWHFIT so that the convolutlon-by-FFT can be used as described 
later. For our recent work, the derivatives In FILLO are computed by making a step- 
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change In each of the estimated fit parameters and recalculating the fit function In 
FILLV. and the array SVPRM In SWHFIT contains the step-sizes currently used. 

SUHFIT has been a general research tool for a variety of waveform work and we have 
found It useful to be able to change the order and number of parameters being fitted 
without having to change the program source code; an array JOROR allows this as described 
further by comments In the Appendix 6 source listing. 

It has sometimes been useful to limit the size of the corrections L allowed In equa- 
tion (14) for the first several Iterations to prevent the Iterative! process from junking 
completely outside the region of convergence In the event of a particularly poor first 
guess for the parameters G^. Subroutine SWHFIT allows 1/5 the correction for the first 
Iteration, 2/5 for the next, and so for^h until the fifth and subsequent Iterations when 
the full correction Is made. 

Another feature of SWHFIT Is Its provision for parameter constraints as suggested to 
us by Dr. William Wells of Washington Analytical Services Center, Inc. (of EG&G, Inc.). 
These constraints are added to the diagonal elements In the matrix A, with each constraint 
being approximately the Inverse of the appropriate fit parameter's Input variance estimate 
In cases In which some a priori Information Is available on the variation expected In the 
parameters to be fitted. 

A general simulation program, to be described only briefly later In this report, has 
been used as an “experimental" means of setting the constraints used for the SE4SAT 
analysis. The ^hge recipe for limiting A In the first several Iterations was the result 
of earlier data analysis and simulation. The exit criterion, based on the fractional 
change In the sum of residuals squared, has the value 0.001 In the SEASAT work, and the 
maximum number of Iterations allowed Is 30. SWHFIT contains feat’ires which are f.pecific 
to the FILLV function being a radar altimeter mean return waveform. One such feature Is a 
test which bars the possibility of a negative rlsetlme or a negative significant wave- 
height. Notice also In SWHFIT that when the antenna pointing angle Is being fitted, there 
are limits on the range of values allowed. Also In the case that the sum of residuals 
squared Increases after several Iteratlw^'s, the parameter values returned by SWHFIT will 
be those for which the sum of residuals squared had Its minimum. 

For the general Iterative non-linear least-squares procedure described, there Is no 
fundamental guarantee that the sum of residuals squared which Is reached In ; fit Is the 
global minimum and not Just a local minimum. We can only point out that for a number of 
cases tried In the simulation and verification procedure, both with and without noise 
added to the simulated waveform samples, the correct parameters were recovered In the 
fitting procedure. Furthermore, the results seened Insensitive to a range of first-guess 
values of the parameters. 
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MODEL WAVEFORM FOR RADAR ALTIMETER DATA 

In this section I will briefly review the fundamental assumptions underlying the 
model waveform used for fitting to the radar altimeter waveform sampler data. The three 
separate terms to be convolved In the model are then separately discussed. One of these 
terms has, for SEASAT, been impossible to represent by a simple analytical expression with 
the consequence that at least one numerical convolution step Is required In generating a 
model waveform for any particular set of waveform parameters. To reduce the computation 
time required, this numerical convolution step can be performed through Fast Fourier 
Transform (FFT) techniques, and this is further described In a later subsection of this 
section on the model waveform. 


Convolutional Model for Rough Surface Radar Returns 

In 19b/ Moore and Williams (ref. 18) demonstrated that the average radar power return 
from a rough surface and for near-normal incidence could be expressed as a convolution of 
the "transmitted pulse shape" and a term which Included the effects of antenna pattern and 
off-nadir pointing angle, surface properties, and distance. More recently. Brown (ref. 4) 
discussed this convolutional model and, for the assumptions common to satellite radar 
altimetry systems, produced a simplified closed-form expression for the average rough 
surface Inqsulse response function. 

My own work is based on the material In Brown's paper (ref. 4) as a starting point 
and I will assume that the average waveform can be viewed in an "elapsed- time" domain in 
which the actual hardware waveform sampling gates can be labelled in successive time 
units. For instance, the 63 SEASAT waveform samplers are spaced over a time span of 
almost IBS nanoseconds. I will assume that I can use the concept of the "effective" pulse 
and ignore the actual details of any pulse compression scheme actually used. SEASAT for 
instance had an effective pulse width of some several nanoseconds, even though the actual 
signal transmitted from the altimeter was about 3.2 microseconds in duration. 

Since the waveform model is viewed in the altimeter's elapsed- time domain, any sur- 
face elevations must be converted to elapsed times (or ranging times) through the radar 
altimeter's two-day ranging time-to-di stance factor of (-c/2) where c is uie speed of 
light. The negative sign is included because an increase in surface elevation (the eleva- 
tion is positive for directions upward out of the sea surface) corresponds to a decrease 
in the ranging time. 

The following description of fundamental assumptions in the convolutional model for 
rough surface scatter is quoted directly from Brown's paper (ref. 4): 
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'• 1. The scattering surface may be considered to comprise a sufficiently large number 
of random Independent scattering elements. 

2. The surface height statistics are assumed to be constant over the total area 
Illuminated by the radar during construction of the mean return. 

3. The scattering Is a scalar process with no polarization effects and Is frequency 
1 ndependent. 

4. The variation of the scattering process with angle of Incidence (relative to the 
normal to the mean surface) Is only dependent upon the backscattering cross 
section per unit scattering area. a**, and the antenna pattern. 

5. The total Doppler frequency spread (4V^/X) due to a radial velocity between 
the radar and aqy scattering element on the Illuminated surface Is small rela- 
tive to the frequency spread of the envelope of the transmitted pulse (2/T, 
where T Is the width of the transmitted pulse)." 

Over the ocean surface, all of the above assumptions are generally satisfied. 

However, we must always be careful In selecting the averaging time to Insure that surface 
statistical homogeneity Is satisfied. For land scatter, the situation Is somewhat differ- 
ent. and some of the above assumptions may be violated. 

Under these various assumptions above, the model mean return waveform W(t) 1s given 
by the convolution of three terms. 

W(t) = PpjCt) * q^(t) * s^(t) . (15) 

Ppg(t) is the "flat-sea" impulse response function which Includes the radar antenna beam- 
width and pointing angle effects. q^Ct) is the radar-observed surface elevation proba- 
bility density function, and s^.(t) is the radar altimeter's point-target response func- 
tioi.. The convolution of these three terms is represented schematically In Figure 3. and 
each term will be described in following subsections in this report. 

For the SEASAT altimeter the waveform samplers actually measure signal plus noise 
(ref. 9), so an average noise baseline parameter, b. Is needed In a model waveform for 
fitting to the sampled data. Also, an overall amplitude parameter. A, is needed. Finally 
the altimeter's altitude tracker moves the entire set of waveform samplers back and forth 
In time relative to the true position of the mean return waveform, so a tracker-related 
term T(t) must be Included to describe the tracker time-jitter probability density func- 
tion as well as a possible track-point shift, t^. The result of these additional terms is 
that the final model waveform to be fitted to averages from SEASAT waveform sampler data 
Is given by 


W(t) » A Ppg{t) *qg(t) * T(t) + b . 


(16) 
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Figure 3. General radar altimeter mean return waveform as convolution of 3 terms. 


Ue have assumed In our SEASAT work to date that the tracker jitter Is small enough In 
width compared to the q (t) that we could Ignore any part of T(t) except for the possible 
track-point shift, t^. In this case T(t) takes the form of a shifted delta function, 

T{t) - 6(t-t^) . 07) 

We may find on detailed Investigation of the tracker Jitter about the true waveform posi- 
tion, that the delta function In T(t) above may need to be replaced by some other form 
such as a skewed Gaussian whose width and skewness are both functions of the significant 
wavehelght and of the AGC level, but this Investigation has not yet been finished. 

The final model waveform, from equation (16) above. Includes six parameters which are 
varied In least-squares fitting the Input waveform sampler data. These six parameters are 
listed In the table below. Figure 1 In the Introduction to this report showed qualita- 
tively how the model waveform changes as a result of changes In these six parameters. 
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In the special case that the s^(t) and q^(t) In equation (16) can each be represented 
by a skewed Gaussian form and that the T(t) has the simple shifted-delta-functlon form of 
equation (17), I have described In ref. 5 a severa1«term analytical expression for W(t); 
an earlier report of mine (ref. 6) was also based largely on this analytical result. We 
have more recently decided that, given the SEASAT s^(t), there Is no way to avoid a 
numerical convolution whenever equation (16) Is to be evaluated for any set of values for 
the parameters In Table II, and we use FFT techniques to speed the numerical convolutions 
In our computer. 

TABLE II. PARAMETERS VARIED IN FITTING MODEL WAVEFORM TO SAMPLER DATA. 


Parameter 

Symbol 

Parameter Name or 
Description 

Term Within Eq. (16) In 
Which Parameter Appears 

A 

Amplitude 

— 

^0 

Time-Origin (or track-point) 

T(t) 

°s 

Sea Surface rms Elevation 


b 

Average Noise Baseline 

— 


Sea Surface Skewness 

qj(t) 


Attitude Angle (or off -nadir 
angle) 

Pps(^) 


"Flaf'-Surface Average Impulse 

Response Function 


This function describes the average power return for a delta-pulse scattered off a 
"flat" ocean; the quotation marks on flat are to remind us that the sea must be rough at 
the centimeter scale for the Incoherent, rough-surface scatter theory to apply, but that 
any surface roughness effects at tens of centimeters or greater (i.e., at the oceanog- 
rapher's level of interest) are described elsewhere in the qj(t) term. From (ref. 4), the 
flat-surface average Impulse response ^unction is given by 

Pfstt) » exp(-6t) iQ(t^B) U(t) , (18) 

in which 6 is a constant (not the 6-function), 

6 » (4/y) (c/h) cos(2fJ ( 19 ) 
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and 


e - (4/r (c/h)* s1n(2E) . (10) 

with U(t) being the unit step function) and l.(t^) being a nodlfled Bessel function. In 
equations (19) and (20) c 1$ the speed of light* h Is the spacecraft altitude* C Is the 
absolute cff-nadir pointing angle (or attitude angle)* and y Is an antenna beaimrldth 
parameter defined as In Brown's equation (ref. 4) by assuming a Gaussian approxImaMon to 
the antenna gain for an angle e off the axis* 

6(e) » 6 ^^ exp[-(2/Y)s1n^e] . (21) 


If 6^ Is the usual antenna beamwldth* the angular full-width at 1/2-power points* then 


4 

Y 


in 4 

s1n^(ew/2) 


( 22 ) 


Notice that the time In equation (18) Is zero at the Instant of the first non-zero power 
return to the altimeter. 

The amplitude term In equation (18) above contains several other constants: 

c o®(0) e 2 

K = T — 5 e*p(- r Sln‘0 , (23) 

® 4(4.)^Lph3 ^ 

where 

Is the radar wavelength* 

o** (0) Is the ocean surface backscattering cross section at normal Incidence* 

Gq Is the radar antenna boresight gain* and 

Lp Is the two-way propagation loss over and above the free-space loss. 


In the SEASAT altimeter* the radar return signals are normalized by the automatic gain 
control (AGC) circuit and we ignore all the individual terms within A^ in equation (23). 
Instead A^ is treated as a simple amplitude scaling constant which Is set equal to one for 
convenience since the A constant in equation (18) will provide the necessary amplitude 
variation. For future work* note that the combination of fitted A and C plus the altim- 
eter's AGC data would lead to a revised estimate. 

Let's emphasize here that the flat surface inq^ulse response function given by equa- 
tion (18) above is an average-power quantity over maqy pulses and has no meaning on an 
Individual radar return waveform. Indeed, as Moore and Williams emphasize (ref. 18)* the 
range from a signal level exceeded 5% of the time to a level exceeded 9S¥ >of the time Is 
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about 18 dB for the underlying Rayleigh distribution In amplitudes. Arvy Individual pulse 
Is likely to look very different from the mean described In equation (18), 

A sample subroutine 6TFSR for evaluating the Pp^Ct) of equation (18) Is provided In 
Appendix D. In evaluating the 1^ In (18), polynomial approximations from Abramowitz and 
Stegun (ref. 19) are used. In the following t » z/3.75. 

If -3.75 1 2 <. 3.75, find Iq(z) from 

Iq(z) » 1 + 3.5156229t^ + 3.0899424t^ + 1.2067492t® 

+ .2659732t® + .0360768t^° + ,0045813t^^ . (24) 

If 3.75 < z < «, find Iq(z) from 

[zW(-z)] I^j(z) = .39894228 + .01328592t*^ + .00225319t‘^ 

-.00157565t“^ + .00916281t"^ - .02057706t‘^ 

+ .02635537t‘^ - .01647633t"^ + .00392377t“® . (25) 

For digital computer use, equations (24) and (25) have to be rewritten In nested form to 
avoid round-off errors from the higher powers in t, and subroutine GTFSR in Appendix D 
contains the rewritten I^ po1ynoi<i1a1 approximations. 

Finally, the possibility of a "negative angle" has been built Into the flat-sea 
Impulse response function Ppjvt) In subroutine GTFSR of Appendix 0. The angle C Is 
supposed to be a magnitude only and thus a negative angle has no possible physical 
meaning. Notice that the fastest plateau decay possible In equation (18) Is when C ” 0 
in which case 6 = 0, 6 * {4/y)(c/h), and the Pp 5 (t) becomes 

Pps(t) = AjjU(t) exp(-6t) . (26) 

There are, however, experimental situations In which an even more rapid plateau decay Is 
needed. There are at least two such situations: a) the actual ^ Is very nearly zero and 
the noise character of the waveform sampler data will lead to plateau decay rates varying 
both plus and minus around the decay rate In equation (26) above, or b) erroneous waveform 
sampler gain calibration data are used so that the decay of the plateau Is apparently more 
rapid than allowed In equation (26). 

For this reason, a branch Is Included In subroutine GTSEA which evaluates equation 
(18) normally for K > 0 but which evaluates equation (26) with a changing "effective" 

6 > If C < 0. The magnitude of the change In 6 is related to the (negative) magnitude of 
C, and the scaling constant on the 6-change Is chosen so that the average of a set of 
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simulated C ” 0 cases Is zero, «Hth Individual C estimates In this simulation uork varying 
positive and negative. 

A negative ^ result from fitting real Meveform sampler data Is a flag or a warning 
that there could be waveform modeling errors (or calibration errors); If the C estimates 
are only sometimes negative and are consistent with the expected standard deviation of the 
C-estImatlon process, then the negative signs pose no problems and simply make the aver* 
ages of sets of C>est1mates come out closer to zero. The expected standard deviation of 
the C-estlmates can be gotten from a simulation program as sketched later In this report. 


Radar*Observed Sea Surface Elevation Probability Density Function 

The radar*observed surface elevation probability density function q^Ct) describes the 

vertical extent and distribution of the ocean surface radar scattering process, and the 

elevation* to- time transformation Is made through the (*c/2) conversion factor. The phrase 

"radar-observed p.d.f." Is used as a reminder that the effective p.d.f. may possibly 

differ from the geometric surface elevation p.d.f. By geometric surface, I mean the sea 

surface as seen by optical Instruments. Our various waveform analysis work (refs. 6, 7, 

and 8) has assumed the one specific p.d.f. form given In the next paragraph, but the 

reader should be aware that there are some alternative p.d.f. forms being proposed. 

Jackson (ref. 20) proposed one modified s^(t) form based on theoretical work assuming 

Infinitely long-crested waves, while Lipa and Barrick (ref. 1) have proposed a somewhat 

2 

different form which Includes a "height-slope cross-skewness" parameter. My approach has 
been to proceed with one particular often-assumed surface elevation p.d.f. with the 
possibility of later substituting a different form at such time as the other research on 
this problem seems sufficiently convincing. Subroutine 6TSEA In Appendix 0 contains the 
p.d.f. described In the next paragraphs, and the possible future use of different func- 
tional forms can be accomplished simply by suitable rewriting of this one subroutine. 

The radar-observed surface elevation p.d.f. q^(t) used here Is assumed to be the 
skewed Gaussian form given In the time domain by 



where Is the surface rms wavehelght, X 1s the surface skewness, and H, Is a Hermite 
polynomial. The Hermite polynomials for general argument z needed In this paper's dis- 
cussion are 


Hjtz) - 





(28) 


(29) 


H^(z) ■ 2^ - 6z^ + 3 , 


and 


Hg(z) = 2® - 152^ + 45z^ - 15 


(30) 


All satellite altimeter measurements of sea surface roughness have assumed that the 
significant wavehelght (SWH) Is exactly four times the rms surface elevation, and since o^ 
Is In time nits, the SWH In distance units Is 


SWH = 4(c/2)o^ . 


(31) 


The X In equation (2^) Is the surface skewness. Skewness Is a non*d1mens1onal 
quantity, the ratio of a third central moment to the 3/2 power of the second central 
moment, so no specific use of the (c/2) distance-time conversation factor Is needed. It 
Is important to realize however that the sign of the surface skewness In the spatial 
domain Is opposite to the sign In the (ranging) time domain. 

The skewed Gaussian form in equation (27) above was used by Pierson and Mehr (ref. 

21) in Skylab discussion and by Walsh (ref. 22) in GEOS-3 data as well as In our SEASAT 
work (refs. 6, 7, and 8); it Is a low order case of a general probability function by 
Longuet-Higgins (ref. 23) for a random variable that is weakly nonlinear. This form Is 
the result of taking the first two terms only in a general Gram-Charlier series (ref. 24). 
Recent work by Huang and Long (ref. 13) proposed that equation (27) above Is not adequate 
for the surface elevation density and that additional terms must be Included In the square 
of the skewness and in the kurtosis, k^, and that the proper surface elevation p.d.f. to 
use Is 


Qs(t) = 


1 + 


r ”3 


(t/a^) 


n “4 


(t/Oj) 




(t/Og) 


X exp[- ^ (t/o^)2] . (32) 

✓?7T Og 

2 

We chose to ignore the X^ and terms, feeling that these were probably less important 

in the ocean than in the laboratory. Also the magnitude was likely to be less than the 

2 

noise in our case. (Figure 6 of Reference 5 shows the effect of the x^ term for example.) 
While our work to date has been based on the simpler q^(t) of equation (27), it is a 
simple matter to change subroutine GTSEA to the equation (32) form If desired. 
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There are two reasons for strong Interest fn the sea surface skewness estimates X^. 
First, this quantity enters some of the theories concerning possible differences between 
“radar-observed” and "actual" surface elevation p.d.f., the so-called electromagnetic bias 
In radar altimeters, as in Jackson's paper (ref. 20) for instance. Second, the Huang and 
Long work (ref. 13} provides the connection between and X^ and two other parameters § 
and X^. The significant slope S and the dominant wavelength X^ are the two Important 
parameters In a new unified two-parameter ocean wave spectral model of Huang £t al^. 

(ref. 14), and thus there Is a possibility of obtaining surface spectrum Information from 
radar altimeter measurements. 


Radar Altimeter Point-Target Response Function 

The radar altimeter system point-target responsv^ function s^(t) Is the system's 
effective transmitted pulse as sampled by the waveform sampler set. This point target 
response Is primarily the effective transmitted pulse shape but also Includes effects of 
the receiver bandwidths. The word "effective" above Is the verbal rug under which Is 
swept the actual details of pulse compression and 1nq>1ementat1on. 

The only Information on the SEASAT s^(t) Is provided by the altimeter Calibration 
Node 1 In which a portion of the transmitted pulse Is fed back through the receiver to the 
waveform sampler set. The SEASAT altitude tracking and SUH estimation In the altimeter 
hardware were based on an assumed pure Gaussian s^(t) with a full-width at 1/2-power 
points of 3.125 nanoseconds. MacArthur (ref. 9) describes a ground-based data processing 
correction to SWH to allow for the ( — ■ ^ ^ )^ shape which the SEASAT s^(t) Ideally had (as a 
result of the particular pulse compression technique Implemented). I have described In an 
earlier report (ref. 6) the actual s^(t) which Is not symnetrlc nor representable by any 
simple analytical expression. Figure 4 from (ref. 6) sketches this SEASAT s^(t), based on 
SEASAT Calibration Mode I Step 9 data, and this figure provides an estimated value every 
1.5625 ns. This Is the spacing between the five sampling gates at the center of the 
waveform sampling gate set, while the spacing between all of the rest of the sampling 
gates Is 3.125 ns. We have assumed zero values at every half-gate position between the 
3.125 ns-separated samplers, based on the assimed nulls In the ideal ( ^ ~ " )^ response. 

In SEASAT the fundamental 1.5625 ns spacing (of the sampled values on all the terms 
to be numerically convolved) Is fixed by the limits on our knowledge of s^(t) and on Its 
densest sampling being at the 1.5625 ns separation at the center of the sampling gate set. 
In subroutine FILLV of Appendix C, the subroutine which evaluates the model waveform, the 
s^(t) Information Is carried In array SYS(514) and the two integers NSYS and NSCTR. NSYS 
Is the point after which all remaining SYS values are zero (clearly NSYS<514), and NSCTR 
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Figure 4. SEASAT radar altimeter system point-target response from 
CAL Step 9 after gain bias corrections. 

is the index of the peak or center at the sampled s^(t). NSCTR is used in a time shifting 
step which in effect sets the peak element within SYS at the time zero. 

If one can assume that the s^(t) has a skewed Gaussian form like equations (27) or 
(32) of the proceeding subsection, then my Reference 5 provides a several-term analytic 
expression for the model waveform and subroutine FILLV of Appendix C could be replaced by 
a coded version of that expression. This may be useful for simulations of future altim- 
eters, but I expect that any actual altimeter will have a s^(t) at least as complicated as 
the SEASAT case in Figure 4, and that the full convolution procedure of this present 
report will again be necessary. 
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Numerlcfll Convolution by the FR Method 

The preceding discussion In this section on model waveforms has emphasized the need 
to carry out a numerical convolution because of the actual SEASAT s^(t) form. Moreover, 
lay discussion of least implied the desirability of being able to change the 

functional form used for this surface elevation p.d.f. This In turn means at least two 
numerical convolutions for a set of model waveform values. Then the step*est1mat1on of 
derivatives within the fitting process will require four more model waveform sets to be 
generated If six parameters are fitted (In general, six parameters being fitted would re- 
quire six more model waveform sets to be generated, but In this problem we know already 
the derivatives of equation (16) with respect to the parameters A and b). Finally the 
general fitting program will require from five to ten Iterations through the waveform- 
pi us-derIvative loop. In brief, there are many convolution steps required to estimate 
waveform parameters for one set of averaged waveform sampler data, and ways to speed this 
processing must be found. 

One way to speed the convolution Is to take advantage of the convolution theorem In 
Fourier transform theory. This theorem asserts that If one has two time functions p(t) 
and q(t) whose Fourier transforms P(f) and Q(f) exist, the Fourier transform of the convo- 
lution p(t)*q(t) Is the product P(f)Q(f). Consequently, the convolution p(t)*q(t) can be 
performed by: a) taking the Fourier transforms of p(t) and q(t) to obtain P(f) and Q(f); 
b) forming the (complex) product [P(f)Q(f)]; and c) taking the Inverse Fourier transform 
of [P(f)Q(f)l to obtain the desired answer. This procedure also can be used for two sam- 
pled time functions using the discrete Fourier transform (DR) and the Inverse DR, pro- 
viding that the two functions can be set to zero outside a region of Interest. The compu- 
tational efficiencies of the Fast Fourier Transform (FR) result in the above-sketched 
convolution through transforms technique being faster in many cases than the simpler 
straightforward numeric/ 1 convolution, and we have found by trying both methods that there 
Is a speed advantage to convolution by transforms within the model waveform computation. 

There is a large body of literature concerning Fourier transforms, DRs, FFTs, convo- 
lution, and so forth; the textbook by Brigham (ref. 12) Is a source which is qu' te read- 
able and accessible, and upon which the discussion In these paragraphs is based. The 
particular FR and inverse FR programs used here for model waveform computations were the 
subroutines FFA and FFS and their associated subroutines from a published package (ref. 

25) of digital signal processing programs. FFA is a radix 8-4-2 FR routine for a real 
data sequence which performs as many base 8 iterations as possible and then performs one 
base 4 or base 2 iteration If necessary, and FFS Is the radix 8-4-2 inverse counterpart 
of FFA, Th.» inclusion of the radix 8 steps in FFA and FFS leads to longer Fortran source 
code but sh'/rter conyiutation times than for other simpler FR routines, and additional 
timing discussion Is provided by Reference 25. 
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The FFT-convolutlon will require that the total array length Nj» used In each term to 
be convolved, be the same and In general there Is the FFT requirement that Nj be some 
power of 2, Nj - 2*^. The smaller the value of M, the faster the computation, but M must 
be large enough to avoid overlap effects as discussed In ref. 12, Suppose again that two 


(discrete) sequences p(t) and q(t) are to be convolved and that p(t) has non-zero 
sampled values while q(t) has such values. There will be (Nj-Np) and (Nj-N^) succes- 
sive zoro values In the arrays for p(t) and q(t) respectively. Then the requirement that 


Nj > Np<»^N^-l will avoid the overlap problem. If there are Uiree sequences p(t), q(t), and 
s(t) with numbers Np, and N^, the requirement on Nj becomes Nj > hp**'N^‘^N^-2. Obviously 
we are assuming the same saniple time Intervals L for p(t), q(t), s(t) and for the final 


answer. 


In the specific SEASAT waveform problem, the sampling time Is 1.5629 ns, the 
spacing of the gates In the center of the set of waveform samplers. The total width 
spanned by the sampling gates Is around 120 A^ and the flat sea response function Pp 5 (t) 
must be allowed at least this width of non-zero values so Np 120. To allow for up to 
20 m SWH Implies that upwards of 200 Is needed to represent q^(t) out to the ± 4o^ 
level, so 200. Then the s>;v,em point-target response function has a width upwards of 
100 A^, so Np » 100. This leads to a total Ny — 120-*-200*100 • 400, and the next hinhest 
2^ Is at M > 9 for Nj > 512. The actual widths used vary somewhat from this Illustrative 
example, but Nj • 512 Is the total array length used. 

The model waveform evaluation subroutine FILLV In Appendix C uses 512-point arrays 
for the Fp^(t), q^(t) and s^(t) and 512-po1nt FFT evaluation (the actual array dimensions 
are 514 because of the details of the FFA and FFS subroutines used). The time-shift term 
T(t) of equation (17) Is not written as a separate term In the convolution, but Is accom- 
plished through the Fourier transform time-shift theorem (ref. 12); a shift In the time 
domain leads to a phase shift In the frequency domain, so that for example If p(t) and 
P(f) are Fourier transform pairs, then the Fourier transform of P(t-t^) Is P(f) exp (-j 
Ziift^). A similar property applies to the discrete Fourier transform, and FILLV takes 
advantage of this and Includes the phase shift In the complex product In the transform 
domain before doing the Inverse FFT to bring the model waveform back to the time domain. 
Details of this are obvious In the FILLV source code In Appendix C. 

Finally, the convolution via FFT method described above was tried Initially as an 
experiment, to see If the computation speed was Improved over the straightforward convolu- 
tloi of two time functions; there was an Improvement of nearly an order of magnitude, so 
the FFT method was kept In our waveform fitting. Adding one more convolution for a total 
of three terms convolved then became simple and cost less than twice the cime for two 
terms only to be convolved via FFTs, and this allowed us to separate the problem so that 
P-c(t), q.(t), and s„(t) were entered by different subroutines. This modular structure 
will easily permit future possible changes in the functional form used for q^(t). Also as 


23 


m note later, most of the time In our wavefom fitting ceeiMtetlont It talwn by the 
computer doing FFTi, and this makes attractive the process of using one of the available 
array processors In any future computer system which spends an appreciable part of Its 
time doing waveform fitting for modeled parameter recovery* 


AOOITIONM. PROCESSING CONSIDERATIONS 


Use of Calibration Data 

If the set of SEASAT waveform sampler measurements Is 1 ■ 1 to 63, It has been 
assumed that there Is a set of gain corrections g^ and offset (or shift) corrections s^ 
such that a corrected set of sampler values Y^ Is obtained from 

Y^ ■ g^(y^-s^) , 1 ■ 1 to 63 . (33) 

The offsets s^ were to be obtained from the altimeter's standby mode, under the assumption 
that they were generated somewhere In the samp 1er» to- telemetry system Interface. The 
gains were to be obtained from the altimeter's calibration mode II In which a uniform 
signal level Is presented to the Input of the waveform sampler set so that the output dif- 
ferences can be Interpreted as gain variations from which a set of gain correction factors 
can be obtained. We have already described (ref. 6) some of tite problems with the gains 
g^ derived from Calibrate Mode II of the altimeter and, more recently, we found that there 
was an apparent gain change correlated with whether the altimeter's AGC word Is above or 
below 32 dB. The SEASAT altimeter's Calibration Node I supplies the sampled point- target 
response function s^(t) but even here the sampler gain calibration problem cause some 
concern. 

There 1$ work yet to finish on the SEASAT smnpler gain problems, and no final values 
for g^ and s^ In equation (33) can be given here. For purposes of this report. It will be 
simply assumed that any necessary waveform sampler calibrations have already been applied 
to the averaged smple data prior to the •'adar altimeter waveform modeled parameter 
recovery procedure. 


Separate Subroutines for Terms Within Waveform Convolution 

I want to emphasize here the <tes1rab111ty of keeping separate waveform terms In 
separate subroutines because of the freedom and flexibility In then changing any one of 
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thtse terms. Once there Is a numerical convolution step required In a model waveform (and 
I cannot Imagine that this will change In future altimeter systems)* then the use of an 
FFT-convolutlon process means that adding an additional convolution will not even double 
the processing time. For example, Reference 5 gives the answer for the convolution of the 
flat'Sea Impulse response function Fp^(t) and the surface elevation density function q^(t) 
for the particular forms of Pp 3 (t) and q^(t) which I am now using. Yet I prefer to do 
this particular convolution numerically too (through FFT again) In order to preserve for 
future work the ability to try different functional forms for 9^(t), forms which would 
Invalidate the results from Reference 5. 


Computer Running Time Conments 

Table III below lists typical computation time taken by the general six-parameter 
estimating Iterative least-squares waveform fitting program on the Data General ECLIPSE 
S-200 computer system In Building E-106 at Wallops Flight Center for three different 
methods of generating the model waveform and Its derivatives. 


TABLE III. 

COMPARISON OF TYPICAL SEASAT WAVEFORM FIT TIMES ON S-200 COMPUTER 

Method 

Waveform Generation 

Derivative Generation 

Waveform Fit Times 

“Full 

Analytical" 

Analytical, from 
Ref. [5] 

Analytical derivatives 
of Ref. [5] expressions 

1-2 seconds 

“Analytical 

Numerical" 

[Pps(t)*qs(t)3 
Analytical; then 
Simple Numerical 
Convolution with s^(t) 

Step-Change In Each 

Parameter, Reevaluate, 

Fa Waveform 1 
” [a Parameter] 

2-3 minutes 

"Full 
Numerical 
by FFT" 

FFT for Pp 5 (t).q^(t) 
s^(t). Then Complex 
Product and Inverse 
FFT 

Step-Change ... 

15-25 seconds 


In planning future altimeter data processing systems, the 20 seconds per fit of the 
“Full Numerical by FFT** entry in the table above would seem to pose a problem since the 
waveform data will probably be averaged and waveform fits done for typically every 2.5 
seconds of data. The factor of 8 between processing time and real-time would create an 
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Intolerable data processing bottleneck, but the existence of add-on hardware array proces- 
sors (and their associated software) suggests a way out of the problem. In Table III, the 
additional time of the ‘'Full Numerical by FR“ compared to the "Full Analytical'* Is pri- 
marily spent doing FFTs In the Fortran programs. Let's assume that all the rest of the 
waveform fitting processing time could be gotten below one second. This Is probably a 
reasonable goal given the relative clumsiness of the Fortran source code written for the 
waveform fitting programs which were developed as research tools and not as production- 
oriented programs. Then let's estimate how many FFT and Inverse FFT operations there are 
In one typical waveform fit, and estimate the operating time required by one representa- 
tive array processor (ref. 26) which can be interfaced to Data General computers such as 
the ZCLIPSE S-200 on ^iiich Table III is based. 

To evaluate a moc-'' ••aveform, subroutine FILLV In Appendix C does two FFTs (on Pp 3 (t) 
and q (t), with transfori: of s (t) already having been performed once in the program 

^ f 

initialization), then three complex multiplications (multiplication of the transforms of 
PF$(t), <l 3 (t), and s^(t) and of the T(t)-generated phase shift), and then one Inverse FFT. 
Thus in effect three FFTs and three complex multiplications of 512-point array quantities 
are required to generate one set of model waveform samples. The MSP-2 array processor 
from Computer Design and Applications, Inc. (ref. 26) does a 512-point real FR in 3.26 
milliseconds and a 512-point complex product In 1.28 ms so that the total time per model 
wavefom) is about 3(3.26)-t-3(l .28) ms or about 13.6 ms. Then the derivatives for the six 
parameters fitted will require generation of four more model waveforms for step changes in 
different parameters in subroutine FILLD in Appendix C, so that one pass through the 
Iteration loop in the least-squares fitting program uses a total of five waveform-gener- 
ating steps or 5(13.6) - 68 ms. Then the average number of iterations in the waveform 
fitting process is generally eight or less so that the total waveform generating time in 
one complete waveform fit is around 8(68) - 550 milliseconds. If the array sizes dis- 
cussed were to Increase to 1024 instead of 512 this time would become about 1.16 seconds, 
slightly more than twice 550 milliseconds, based again on the NSP-2 processor literature 
(ref. 26). 

Thus the entire waveform fitting procedure for one set of waveform sampler averages 
could probably be carried out in under two seconds of computer time on an ECLIPSE S-200 
computer with an array processor. I use the Computer Design and Applications MSP-2 
array processor only as an example to show that a 2.5 second average of waveform data 
could be processed in less than 2.5 seconds of computation. Other array processors are 
available for minicomputer or large mainframe computer systems at varying prices, speeds, 
and features. At least one unit available for ECLIPSE computers runs more than twice as 
fast as the MSP-2 processor In the example above. No array set-up or transfer times or 
other system overhead times have been included in the above examples, and the possible use 
of an array processor deserves a more detailed study for future altimeter waveform proc- 
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essing. Final ly^ It should be noted that the computation times required by maqy of the 
current 32*b1t "super-minicomputers'* Is considerably less than the computation time of our 
ECLIPSE S-200 system; It might be possible that a newer computer would operate fast enough 
that no array processor Is needed, but this again Is a further Issue for detailed study. 


Simulation Prograun for Assessing Operation of Fitting Program 

The model waveform generating procedure has been discussed earlier. Verification of 
correct operation of this program Is fairly simple; one generates a waveform for a given 
parameter set and checks that the result matches what the theory predicts. One way to 
test the SEASAT waveform program FILLV of Appendix C Is to write a small test program to 
produce a skewed Gaussian characterized by a rms a and X and to then take 1.5625 ns- 
separated samples of this function to form a point-target function s^(t) In array SYS; 
then the waveform generated by FILLV for these SYS must agree closely with values calcu- 
lated from the method described by Reference 5. This has been done for the SEASAT model 
waveform used here und the verification was satisfactory. 

Another level of checking Is needed as well to verify that the waveform least-squares 
fitting program recovers the correct parameters from the input waveform samples. For this 
check, a simulation program was developed which generates model waveform samples for one 
particular set of waveform parameters and then uses those sample values as input to the 
waveform fitting program; a comparison of the parameters recovered and the Input param- 
eters will then reveal any errors or biases induced by the fitting program. This simula- 
tion program Is also extremely useful In determining the fundamental limits to parameter 
recoverability, limits due to the fundamental noise-like character of the Individual radar 
return waveform. 

Figure 5 sketches the simulation program used here to Investigate waveform parameter 
recoverability. The average of N individual returns, for any particular specified time- 
point on the waveform, is expected to have a Gaussian distribution because of the central 
limit theorem. The Individual return itself is described by an exponential distribution 
whose mean and standard deviation are equal, so the standard deviation of the N-average 
return is expected to have a standard deviation equal to the mean divided by v^, and this 
noise modeling Is built into the "noise loop" shown in Figure 5. The results of consider- 
able experimenting with this simulation for the SEASAT work showed that the waveform fit- 
ting program operated with no appreciable biases and appeared tolerant of a reasonable 
range of first-guess parameter values. No further details will be presented here as this 
is the subject for a future separate report. The important point for this report is that 
this sort of simulation program must be produced and exercised for each different radar 
altimeter's particular conditions, and that work on waveform processing software Is not 
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PROBLEM PARAMETER LOOP 


complete until the simulation program has exercised the waveform modeled parameter re- 
covery programs. 



Set: NR > I Individual returns to be averaged; number of waveform parameters 
to be fitted for; range of waveform samplers to use In fits; NWF > 
number of noise waveforms to select; other Initial conditions; and 
the Initial, final, and delta values for the problem parameters 
such as SWH, pointing angle, and skewness. 


Select next set of problem parameters such as SNH, and 

Compute the 63 expected sampler values Y , 1»1...63, for 
these problem parameters. ”1 

Zero the noise loop sunning locations, set NFIT»0. 

Do llnefit, print recovered waveform parameters for the 
zero-noise limit case . 


NFIT ♦ NFIT+1 

Select a set of rzndom numbers R^, 1<‘1....63, with a 
unit mean and standard deviation 'of l/4iR . 

Form the average-of-NR return samples Y. by Y.=R.Y_ , 

1=1. ..63. ' ' ’ ®1 

Do llnefit, store recovered waveform parameter estimates. 


NFIT < NWF? ^ 


a. 

s 

ui 

(/> 

o 


Yes 


< 


Form mean and standard deviation for each of waveform fit 
parameters, and print out these results 


More Input problem parameter cases to run? 





Figure 5. Simulation program for assessing waveform parameter recoverability. 
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APPENDIX A 


GENERAL DISCUSSION RELATED TO WAVEFORM FITTING SUBROUTINES 


The following four appendixes In this report provide source listings for various 
subroutines In the waveform fitting procedures and for SEASAT-1 data. The main program, 
which Is not provided here, handles obvious details of data Input, averaging, and call* 
bratlon. There are also some Initial set-up details handled In the main program, setting 
the number of parameters to be fitted, lower and upper sampling gate limits If fewer than 
63 waveform sampler results are to be used, setting the array for the system point- target 
response, and so on. The main program also sets the values of arroy SYS which carries the 
sampled system point-target response function s^Ct), and sets the values of an array GUESS 
which carries first parameter guess for each parameter being determined In the Iterative 
fit within SWHFIT. The array also carries the fixed values of the parameters not being 
fitted If fewer than six parameters are being determined. 

Table A-1 lists subroutines used In the waveform fitting and Indicates where their 
source listings are to be found. Table A-2 shows which of the subroutines are called by 
others; obviously subroutine SWHFIT Is called from the main program. In the Data General 
FORTRAN used, “X" In column one of a source line means the line Is to be Ignored unless 
the /X switch Is used In the FORTRAN compile step. Various debug statements In the source 
listings, with X In column one, call subroutine RLOUT and CPOUT which are not supplied. 
RLOUT and CPOUT are trivial array formatting and printout routines for real arrays and 
complex arrays respectively. 

The calling list for SWHFIT Includes the 63 waveform sampler corrected data averages 
In YIN(63), and sampling gate numbers NLO and NUP. Data for gates from NLO to NUP will be 
fitted, and normally NL0»1 and NUP=63. The calling list array WTY(63) carries relative 
weighting factors and, for the generally-used uniform data weighting all the elements of 
WTY are set to unity. Notice the array STPRM(7) which Is set to the step sizes to be 
taken In producing the parameter derivative estimates by FILLD. Actually the first, 
fourth, and seventh values of STPRM have no meaning since the amplitude (first) and base- 
line (fourth) derivatives are found without making a step-change In parameters one and 
four, and the seventh value of STPRM was Intended for the kurtosis which has not (to date) 
been written Into the function In GTSEA. Notice also In the SWHFIT comments the descrip- 
tion of J0RDR(7), an Integer array determining the order In which the parameters are to be 
found . 

Part of the data communication In these subroutines Is by calling lists and part by 
named COMMON areas. Table A-3 lists the contents of named COMMON area SSM4N, and Table 


A-4 Msts the contents of named COMMON area SYSTM. Two additional named comnon areas CON 
and CONI are used within the FFT package only. 

The Individual subroutines have a number of additional convnents provided In their 
source listings. It Is hoped that these comments, the Tables A-1 through A-4, the remarks 
above, and the the main body of this report should provide the reader an adequate ac- 
quaintance with the Wallops waveform modeled parameter recovery procedures. 


TABLE A-I. SUBROUTINES USED IN MAVEEORM MTOELEO PARAMETER RECOVERY PROGRAMS FOR SEASAT-1 DATA 


Routine Nm 

Where Provided 

Purpose and Oescriptlve Remarks 

SUNFIT 

Appendix B 

Function subroutine to fit up>to-6-paraneter model MBveform to 
corrected averaged Maveform sampler data supplied by main calling 
program. 

SYMINV 

Appendix B 

Symnetric matrix Inversion subroutine used by SWHFIT. 

FIUV 

Appendix C 

For given parameter set, generates model Maveform at the 63 Maveform 
sampler locations. Waveform is generated through FFT-numerical 
convolution of three separate terms. 

FILLO 

Appendix C 

Uses FIlLV. changes each parameter by specified step to form estimates 
of derivatives of model Maveform Mith respect to each of the waveform 
fit parameters. 

GTFSR 

Appendix D 

Produces flat-sea impulse response function, one of the terms con- 
volved in FILLV. 

GTSEA 

Appendix 0 

Produces radar-observed surface elevation p.d.f., one of the terms 
convolved in FILLV. 

FFA 

Appendix E 

Does Fourier transforms in radix 2-4-B FFT package used. 

FFS 

Appendix E 

Does inverse Fourier transforms in radix 2-4*8 FFT package used. 

ORD1 . 0R02, 
R2TR. R4SVN, 
R4TR, R8SYN, 
R8TR 

Not Provided 

Additional subroutines in radix 2-4-8 FFT package not provided in 
this report but available from Reference [25] or as described in the 
comnents in FFA and FFS in Appendix E. 

RLOUT, CPOUT 

Not provided 

Routines for printing out real and complex arrays, respectively, and 
used for debugging only. 
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Debug (optional) Routines RLOUT, or CPOUT to Print Values in Real or Complex Arrays Respectively. 


TABLE A-3. DESCRIPTION OF VARIABLES IN NAMED COMMON AREA /SSM4N/ 


Variable 

A{7) 

AEOIT(2,7) 

CALM 

CNSTR(7) 


CORRL(21) 

GUESS(7) 

ITER 

NA 

JOROR(7) 

SERSQ 

T(63) 

XCNST(7) 


Description 

The seven final values of the up-to>seven parameters fitted to 
Input waveform data. NA • number fitted. (The seventh Is kurtosis* 
not yet written Into 6TSEA. so A(7) > 0 always.) 

Provides for possible lower, upper edit limits for the A(7). 

"Calm-Sea Constant" or effective l-o point-target response width 
used In finding SWH with "analytical" fitting. When numerical 
convolution (via FFT) Is used, main program should set CALM ■ 0. 

Constraints on the fitted A(7) which In effect limit possible 
size of parameter changes In the fitting Iterations. CNSTR 
values are set by a priori Information, once, at start of main 
program. 

Provides Information on correlation of fitted parameters - 
probably of doubtful use In this non-Kncir problem. 

Carries first-guess for waveform parameters to start Iterative 
fit procedure, and carries fixed values for parameters not being 
fitted. 

Indicates # iterations needed to fit data upon return from SWHFIT. 

Number of the A(7) being fitted, with maximum NA > 7 (If kurtosis 
1$ added, for example - current maximum NA Is 6). 

Set by main program. Indicates order of parameters being fitted, 
and provides ability to fit different parameter combinations 
without rewriting source code. 

Sum of squares of fit residuals, upon return from SWHFIT. 

Carries the 63 Independent parameter (sampler time) values. In 
SEASAT, T(l) » -92.1875 and T(63) - 92.1875. 

Provided for transporting different constants, but not used for 
anything significant at this time. 


TABLE A-4. DESCRIPTION OF VARIABLES IN NAMED COMMON AREA/SYSTN/ 


Variable 

NSYS 

NSCTR 

SYS(514) 


SYSIO(21) 


Description 

Value of index I after vfhich all SYS(I) are zero. 

Center-of-response index, the index within SYS to which is 
attached the time zero label. 

Contains the sampled point- target response in the elements 
SYS(1) through SYS(NSYS), with all remaining elements through 
SYS(514) zero ^fore the first call to FILLV. After that, SYS 
contains the 5lZ-point FR transform of the input s.(t) and 
must not be changed. 

A comment line carrying Information on how SYS was determined. 
SYSID is only of interest to the main program. 
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can c 
i3r02 c 
arrj c 
0JJ4 c 
B€Oi c 
oy:-si c 
c 

of-ao c 
Cisys c 
OvJiC c 
cm c 
01*1 2 c 
0C13 C 
02.14 C 
0013 C 
0:716 C 
0,'317 C 

rci 3 c 

0019 C 
L920 C 
0321 C 
0022 C 
5023 C 
0324 C 
002b C 
032C C 
0027 C 

ro20 c 

0123 C 
0330 C 
0321 C 

0332 C 

0333 C 

ccz** c 

0035 C 

0036 C 
5037 C 
0038 C 
OJ39 C 
eHji9 C 
£J41 C 
00*2 C 

con c 
c:u c 

0340 C 
0J4O C 
0347 C 

CJ-.S c 

0J’,9 C 
C3i0 C 
CLbl C 
0.332 C 
C3j3 C 
JOJ 1 C 
OO'.b C 

o-;j 6 c 

0*J;7 C 
f'"">3 C 


SOURCE LISTINGS FOR SUBROUTINES SWHFIT AND SYMINV 


DPCtHAYNEtSV.nFITS.FR REVISED 10/30/10, APPROX. 1316 HRS 

G.S.HAVrCE WTC/DAS 09/26/80 

SU3R0UTINE SWHF1T3 IS EXACTLY SA,*« AS S1kWIT2 EXCEPT THAT 
A NEW C0M10N AREA /SVSTM/ HAS BEEN ADDED A4D THE COMIUNICATION 
WITH FILLV I FILED CHANGED SLIGHTLY. THIS IS SO THAT THE 
SYSTEM POINT-T.««GET RESPONSE DATA FROI FILLV IS AVAILABLE 
BACK AT THE MAIN CALLING PROGRAM. 

S.MALL CHANGE ON 10/30/80 CHANGED WHAT HAPPENED FOR A(3) 

NEAR ZERO. 

aiANGE ON 07/16/80 ADDED ARGUMENT IXDdG TO ROUTINE CALLING 
LIST; A VALUE OF 1 FCR IXDBG SIGNALS PRINTOUT OF EACH 
INDIVIDUAL ITERATiai IN THE FIT. 

THIS VERSION FOR THE 512-POINT TRANSFORM. AND COfTXlN AREA 
SSM4N. 

THIS IS REVISION OF SWHSSN (REV. 07/27/79) NOW INCORPCRaTINC 
FFT-CONVOLUTION METHOD OF GENERATING SEASAT-l FUI4CTI0N VALUES 
AND DERIYATIVES. THE SUBROUTINES NEEDED BV SWHFIT INCLUDE: 
FILLV3.FILLD3, FAST . FSST , FR2TR , FR4TR , FR4S /N , FORDl , 
FORD2,GTSEA2,GTFSR2, AND SYMINV. 

NLO £ NUP ARE NOW IN CALLING LIST OF SWHFIT; THZSc DEFir.’E 
LOWER AND THE UPPER LIMITS TO THE r.T.S'OER OF POINTS BEING 
FITTED. IT IS REQUIRED THAT NLO.GF.I, NUP.LE.63. ATJ) THAT 
N.V.LT.(NUP-NLO). 


AS IN SIHSSN, THE INTEGER ARRAY Ov»DR(7) DESCRIBES THE ORDER 
IN WHICH PAR.1METERS ARE TO BE FITIEO. THE FIRST NA (NA.LE.7) 
ELEMENTS OF JOkDR(.) ARE TO OE INDIVIDUAL (AND DIFFCREND VALUES 
Or 0, WriERt 0 MAS THE FOLLCrJING MEANING: 

0 ■ 1 DES1G4UTES AfPLirUOE 

“ 2 ” TIME ORIGIN 

" 3 " RISE TIME 

“ 4 " EASELIME 

• 5 " SKEWNESS 

“ 6 " OFF-NADIR A^IGLE 

“ 7 * KURTOSIE. 


A 7-PARAriTEP NON-LINE.\R LEAST- fC-U/RES FirTIHG rU.'CTlON SUB- 
ROUTINE SPECIALIZED TO CASE OF SEAb.^T SIG.‘IIF:CANT VA/EMEICHT 
DETERrilNATICM. SUBROUri.'E S\»HKIT INCLUDES POSSIhlLITY OF UVuQUAL 
WEIGHTS FDR THE 63 INP'JT Siri DATA POINTS. WEIGHTING BY INVERSE 
OF VARIANCE ESTIMATES. FOR EXAMPLE. PROCUCES IN EFFECT A 
maximum-likelihood fit. THE INOIOVOR IfR IS 1 ON THE RETURN 
FROM A SUCCESSFUL FIT; OTHEPWISc ICR CAN INOICLYc THE TYPE OF 
ERROR AS INDICATED IN CO^•.'lENTS BELCW. 

THERE'S ALSO PROVISION FOR CALCULATING t0.n£L,V!!0N COEFFICIENTS, 
Kuj ABILITY TO ADJ A PRICkI FIT CC'JiTRAINTS, RESOLTINS FROM 
COriVERSATlCNS WITH BILL WELLS OF VA2CI . INC. 



L’TO'd 

OSO0 

OCSl 

DJC? 

0CIS3 

W64 

f'wes 

(/0S6 

C0S7 

0;i33 

0059 

0370 

0971 

0372 

0073 

0374 

0375 
C'';S 
00/7 
03/C 
037 D 

0330 

0331 
9Ci'.2 
0333 
033 « 
0335 
OC36 

00 J 7 

oj;o 

0.VS9 

0J90 

3091 

0032 

03^3 

0394 

0333 

CSOo 

crs7 

0-393 

C399 

DIOS 

OlfJl 

0,J,2 

DVji 

ai'Vi 

C1J5 

0176 

0107 

o:i;8 

r.-t 

01 10 
OUl 
0112 
«'P 1 
0114 
01 la 
0116 
M! 17 


c 

C SUSROUTINE SWHFIT USES SVK'IETRIC K/sTKlX INVCRSIUN ROUTINE SM11NV 

C TO INVEST SVr.'«TRIC MATRIX XMAT IM PROBLEM. FUNCTICNAL FORM TO 

C EE FITTED IS SUPPLIED BY SUBROUTINE FILLV AND THE NEEDED 

C DERIVATIVES BV SUBROUTINE FILLO. 

C 

F'JNCTION SWHFIT(NL0,NUP,VlN,tfrV,lXD30,lER) 

COMPILER STATIC 
C 

REAL STPRM(7)/1 .0, 0.2S, .20, 1.0, 0.03, 0.03, 0.03/ 

C 

DIMENSION VIN(63),WTV(63),yMAT<7,7),BCOLM(7),PVECr(7), 

1 WT(63),PSYM(7),OSYM(7),MSVM(7),CN8TI(7),AXEEP(7),XKEEP(7,7), 

2 VAL(fi3>,DRV(63,7) 

C 

C 

COMMON /SVSTM/ SVS(S14),N3VS,NSCTR 
C 

COMMON /SSM4.N/ A(7).XCNST(7),T(O),NA.ITER,SERSQ,CALM,C0RRL(21), 
1 GUESS(7>,CNSTR(7).JORDR(7),AEDIT(2,7) 

C 

X WRITE FREE (12) •<15> AT SWHFIT3 IllPUT, WTY VALUES ARE:* 

X CALL RLOUT{12,WTY,63) 

C 

C... SET INITI.U VALUES. LIMITS 

C 

C 

IF (NLO.LT.l) NLO-1 
IF (NUP.GT.63) NUP»63 
C 
C 

c 

ERLIM*0.031 

LXMIT-OO 

1TER»0 

SWMFIT»-3SB8. 

IER--10 

WT1»0. 

SaMIN»l.E36 

sersq=scm:n 

c 

C (USED TO REMOVE BIASES HERE.) NOW CHECK THAT IN?UT VALUES ARE 

C... WITHIN A1L0*'£D LIMITS (RETURN UITH NO FJSTME? SWH LVRK IF 
C... NOT WITHIN LIMITS). AND SUM INPUT WEIGHTINC FACTORS FOR 
C... NORMALIZATION. 

C 

DO 5 I-.N'.O.NJP 
WTI»WTI*UTV(n 
YI»VIN<I) 

IF ((VI.LT.-25.).O«.<YI.CT.300.)) RETURN 
5 CONTINUE 

C 

DO 10 I'NLO.NUP 
19 WT(1)«V/TY(I)/Vai 

C 

c 

c 

C... SET INITI/X FIT PARAMETER ESTIttATES. SET ALL TERMS FROM GUESS<7) 

C.. WHICH IS A£SL-i"£0 TO SET VALUES NOT BEING FITTED .AS WEIL AS 
C.. THr INITIAL GUESSES FC.T THOSE BEING FIi lED. 
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0124 
0126 
0120 
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0130 
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0102 

0133 

0134 

0135 

0131 
0137 
0*35 
0139 

01 ’.g 

B141 

C142 

0143 

0144 

0119 

0145 
P!;"? 
C143 
0149 
0:50 

0151 

0152 
015J 

0154 

0155 
015C 

01-j/ 

0155 
019‘i 
P'6‘9 
0151 
PU2 
P143 
0U4 
0165 
0i6e 
01 CV 
016C 

0169 

0170 

0171 
6172 
C173 
0174 
PI 75 
0176 
017/ 
I.' ?l 


C 

DO 9 I«l,7 

9 A(I)«CUESS(I) 

C 

C— AT THIS POINT IN PROGRAM. SET THE INVERSE-VAJIXANCE CONSrRMNTS 
C-- TO OE ADDED TO 1HE ON-DIAGONAL TER.1S BELOW AT "DO 92* Si‘Alu-‘icNT. 

C— PUT TEST TO AVOID AN INPUT ST/tTlARO DEVIATION ESI I MATE LESS THNl 
C-- /30V: 0.033 
C 

I3K>0 

UO 12 I1>1,NA 
l-OCRDR(II) 

IF (1.EQ.6) IJK-Il 

r... CHECKS WIHETHER POINTING ANGLE IS ONE OF PARAMETERS BEING FITTED, 

C.. SINCE THERE MAY BE CORRECTION SUDSEQUCNTLV TO AVOlO A tlESAIlVE A.'.GLE. 
VI«CNSTR(1)*»2 
IF (VI. LT. 0.001) VI-0.001 
CNSTim-l./Vl 

12 AKEEP(1)-A(I) 

C 

c 

c 

c--- FILLV CtTS THE 63 SA'^PLED VALUES FRO-5 THE FFT-Ca.VOLUTIOH 
C— - PROGRAM. THESE AP.E IN ASCENDING O^OLR IN IIVEPENDcKr VARIABLE. 

C 

6 CALL FILLV(VAL) 

tOLD’O. 

ILR-t 

DO 13 I-NLO.NUP 
VI->VIM(I)-VAL(I) 

13 EOlO-EOLD ♦ VI«V1«WT(I) 

C 

C 

X ITliTE FREE (12) ‘(IS) IN SWSFIT, THE VIN VALUES ARE :* 

X CALL Rl0UTU2.V1N.63) 

X WRITE FREE (12) *<16> AND FIRST-CALL VAL FROM FILLV ARE :* 

X CALL RLOUTU2.VAL.63) 

C 

SCMIN-EOLD 

C 

c 

ir (ixdbs.se. 1) CO to ii 

c 

WRITE FREE (121 *<15> SWHFIT; NA -•.NA.‘ ITCRATICS RESULTS APE:* 
U-RITE (12.14) ITER. COLD. A 

14 FORMATC ITER •‘.IS.’. SERS'J -'.614. 6. *; A( . ) ■ ;*.(* •.7F11.9)) 
C 

C... 5TATL.‘-iENT 11 IC RETURN POINT FOR NEXT ITE.V,TIC.'J IN FITTUk; . . 

C 

11 CONTINUE 

C 

C... ZERO UPPER PART. SVPTCTRIC MATRIX . . 

C 

DO 20 I-l.NA 
DO 15 J-1.1 

19 >3*j4T(O.I)-0. 

20 ECOLFHl)>^0. 

XFRCT-1. 

C— -PUTTIUC A COMXEf.T IN FCLLDWIsr C.*?1D WILL PE7SVE EFFECTS Of 
C-- TMiS AiTtfiPT (2/17/78) AT RED«.>CiKG S!2C OF FiMST FtU *STIPi* 


Li*. 70 
0130 
OlSl 
0132 
013’J 

0134 

0135 

0136 
0107 
UlfIS 

0189 

0190 
0!O1 

0192 

0193 

0194 

0195 
0190 
0197 
0193 

rio«» 

CiZ'3 

C20? 

02<C 

.02J4 

0209 
0200 
0207 

02:2 

0200 

0210 
0211 
tv 12 
021J 
0214 
0'lt 
0210 
0217 
u2)0 

0. M9 
r;220 
0 ??! 
1:222 
0273 
0224 
4>'22u 

02 ?r. 
1 22' 
0223 
0220 

1. ;:o 

0:11 

0202 

02^3 

10:1 

i»22% 

nz’i 

tiv:: 


C— IN PA.V>!ETrR SPACE . . . 

IF (ITER.LT.4) XFRC1 « < 1 . *FLCAT< ITER ) )/6. 

C 

ELIM'-ERLIM'XFRCT 

SEP.SQ-0. 

C 

C— > SUdROUriNE FIUO SETS OP THE (63,7) DERIVATIVE ARRAY BV MAKING 
C— STEPS IN THE VALUES OF THE PARAMETERS A(7>; THE STEP SIZES 
C-— TAKEN ARE CARRIED CV STP?.M(7), AllD THE OR0EP. OF THE DERIVATIVES 
C->- IN 0RV(63.7) IS SET BY 3CRDR. FIUV MUST HAVE BEEN CALLED BEFORE 
C— FILLO; THE MELD VALUE IN VAL(63) FROM THE CALL TO FILLV IS USED IN 
C*— FILLING DRV(63.7). INCIOENTALLV, NUMERICAL OEKIVATIVE IS NOT 
C-— D09E IN CASE OF THE AMPLITUDE AND BASELINE. SO THE VALUES SET IN 
C-— STPrrid) AMD STPRM<4) />RE IRRELEVAT.T. 

C 

C 

C 

CALL FILLD(VAl .DRV.STPRM) 

C 

C_»_ -00 36- LOOP FILLS UPPER HALF OF THE SYMMETRIC MATRIX, ALSO THE 

C-- COLU.-W VECTCR . . . 

DO 33 OP-NLO.NUP 
WTI-WT(JP) 

DV«VIN(J?>-VAL(JP) 

C 

C 

DO 35 0A>1,)UL 
PVCCT(OA)«DRV(OP,OA) 

35 CONTINUE 

C 

DO 30 I«1.NA 
VI«PV£CT(I)*WTI 
DO 40 0*1.1 

40 XMAT(J,I >»XMAT(0,I)*V1*PVECT(0) 

93 bCOLMdl’ACOLMd >*V1»0V 

C 
C 

C--- NOW AOO ON-DIAr.O-.'AL CONSTRAINT ELEMENTS TO THE SYMMETRIC MATRIX; 
C-- THESE APE FROM Tl.E A PRIORI IMVORMATION ON VARIATION EXPECTED IN 
C-- THE Mr.A' ETERS TO BE FITTFD IREF. COiWERSATION WITH BILL WELLS OF 
C-- UASri. ABfyjT C8/2;/78 ). THE CO J STRAIN! IS THE INVERSE OF THE INPUT 
C-- VAP.r/iwE ESTlfUTE. MINIMUM ALLCWEO VARIANCE OF 0.031 
C 

c 

c 

CO 52 1*1. t.A 
OJ«JCRwT(T) 

52 XMAT(I.!)*XMtArd.l)*aiST!(OJ) 

C 

C--- CALL MtRE THE S/.'YtCTRIC MATRIX INVERSION ROUTINE FR«I 0. MCMILLAN 

C 

C 

tX VTITE FREE (12) 'BErOKE MATRIX INYERSIWl. INPUT FUTRIX 

CX DO 779 I*1.KA 

CX 779 V.2ITE (I2.7C0) (XMAT(O.I). 0*1,1) 

CX V83 tC14.6) 

CALL Syii:NV<X,'UT.KA.IFAIL.7.PSVi1,e'Y;4.i1jVM) 

CX IF (IFaIL.NE.O) VV.ITC free d2) '('.ArRIX INVERJItX FAILED . .* 

CX V-IVC ^REE (12) • AFTER F>Tf:IX IN/EAblOM. OUTPUT FUTV.IX 

CX CP 731 I»l,NA 


0239 

CX 7 

91 V.;?1TE (12,783) (X«AT(0,I>, J»l,l) 


C J VJ 

c 



W241 

cy 

WRITE FR.IE (12) ‘SYMIN-V IFAIL SHOyiO ■ 0; IT IS '.IFAIL 

0141 

c 



02 ',3 


IF (IFAIL. NE.0) GO TO 1000 


0244 

c 



0245 

c 



02 4C 


DC 66 1=1, NA 


024? 


:i=JORDa(I) 


0249 


AT' 


0249 


r . 2 0=1, NA 


0233 


* J.LT.l) CO TO 60 


0251 


A'. L M= ACL - l*XMAT (1,0) “BCOLM ( 0 > 


0252 


CO TO 62 


0233 

63 

ACL.1=ACLH*XMAT(0, 1 )»BCOLM(0) 


0254 

62 

CONTINUE 


0235 

C 



0255 


IF (II.NE.6) GC TO 64 


0257 

C 



0253 

C... 

THE 11=6 PARAMETER »S POINTING ANCLE; FOLLOWING TREATMENT IS AD 

0259 

C.. 1 

HOC AND SPECIFIC TO SEASAT-1 CASE (07/25/80). 


0250 


IF (ABS(A(6)).LT. 0.025) ACLM=ACLM/5. 


02-1 


A(6) = A(6) ♦ XFRCT'ACLM 


0252 


IF ( (A(6).LT.-4.).0R.(A(6).GT.4.) ) A(C) » A(6)/10. 


0253 


CO TO 65 


0264 

c 



0263 

64 

A(II) = A(II) ♦ XFRCT*ACLM 


026 S 

C 



?2S7 

C 

. 


0253 

SS 

co7rriNUE 



C 



3273 

C- — 

PUT FOLLOWING STEP IN TC AVOID NEGATIVE RIScTIME. 


9271 

c 



02 72 


IF (A(3'.LT.l.E-06) A(3) = l.E-06 


0273 

c 



6274 

70 

CALL FILLV(VAL) 


3275 

C--- 

RECALCULATE VALUES OF THE SAMPLED WAVEFORM FUNCTION FOR 

THE NEW, 

0276 

C— - 

UPDATED ESTI7UTES Or THE PARA.METERS A(7). 


3277 


SEPSQ=0. 


0278 


DO 175 I=NLO.NUP 


0279 


VI=YirJ(I) - VAL(I) 


3233 

175 

SERSQ=SE?.SQ*VI*yi«WT(I) 


0281 


ITER=1TER-1 


0282 

C 



3133 

C 



3284 


IF (IXDOG.EO.l) WRITE (12,14) ITER,3EkSQ,A 


0285 

c 



0236 

c 



0287 

c 

.. CHECK THAT WE KEEP COEFFICIENTS PRODUCING MINIMUM SUM 

ER.RORS** 

02SC 


IF (SER3Q.GE.SCMIN) GO TO 185 


0289 


DO 180 1=1, NA 


0293 


DO 179 0=1.1 


0201 

179 

XKEEP(0,; )=XrAT(0,I) 


0292 


!I=JCRDR(I ) 


0293 

1C3 

AKECP(ll)=A(II) 


0294 


SOMIN=S;RSO 


0295 

135 

CONTINUE 


0295 

C 



0297 

C 



0' «■: 

C 




B5 


02S3 

0W3 

OZOl 

0302 
C3J3 
03J4 
0305 
IJ306 
0307 

0303 

0309 

0310 

0311 

0312 

0313 

0314 

0315 
0315 

0317 

0318 

0319 

0320 

0321 

0322 

0323 
U324 
0325 
C320 
0327 
0323 
03i'9 
03 30 
0331 
J332 
0333 
033 1 
0335 
03^6 
0337 
0333 
0?3J 
03 -'0 
cm 

0342 

0343 

0344 
03 15 
C3 lb 
03 ,7 
0318 
0349 
0330 
0351 
3352 
033 3 

0354 

0355 
0355 
036 7 

I > fj '3 


C TEST; AVOID TRVING FOR ABS03OLV SMALL RZSIDUAL5 A30UT FIT . . . 

IF (SERSD.lt. 4. 4E-05) GO TO 3CC0 
C 

IF (ITER.GE.LIMI n CO TO 2033 
IF (StRSO.GT.EOLD) GO TO 5030 
C 

IF ( ((EOLD-SEkSQ)/EOLD).LE.ELlM ) GO TO 3000 
190 EOLD=SERSQ 
GO TO 11 
C 
C 

C— FOR THE SEASAT WAVEFORM C/.SE, CHECK IF ERRCRS*<*2 INCREASED; DON’T 
C-— MAKE ERROR EXIT IF (FRACTIONAL) INCREASE IS LESS THAN 10 TIMES 
C— LIMIT. ALSO CHECK ABS. MAGNITUDE OF SERSQ; VALUE LESS THAN .0099 
C-- IS A INDIVIDUAL S£H CATE STANDARD DEVIATION OF ABOUT 0.015 M, 

C— WHICH UE ASSUME TO BE AN ADEQUATE LOWER LIMIT TO THE SEASAT 
C-- SITUATION WHEN USING ONLY LAST 46 GATES. 

5000 IF (SERSQ. LE. 0.0099) GO TO 3000 

IF (ABS(lSERSQ-EOLD)/EOLD).LE.10.»ERLIM) GO TO 190 
C 

1ER»-1 

C---SU1 FRR0RS*«2 INCREASED 
C 

GO TO 2503 
C 
C 

1000 lER=-2 

C— - .’•WTPiX IMVERSICN FAILURE (SINGULAR MATRIX . . .) 

C 

GO TO 2500 
2003 IER=-3 

C— ITERATION COUNT EXCEEDED 
C 

c 

c 

C... IF ITER .CT. 2, FIGURE THAT SttME SORT OF SOLUTION EXISTS. SO SET 
C.. 1ER=1 A'JO RETRIEVE THE MIMMU1- VAI.UE-PROOL'CING SET OF A( . ) AND THE 
C.. RESULTING XMAT(...) VALUES IF NECESSARY . . 

25&iT IF (ITER.LE.2) GO TO 4533 
IER»1 

:r iSERSQ.LT.SOriN) GO TO 3000 
DO 2510 I»1.NA 
DO 2505 J<^1,: 

2505 XMAT(J.I)=XKE£P(J.I) 

II»OORDR(! ) 

2510 A(II)»AnEE=(II) 

StRSG=S&.';IN 

C 

C 

C 

303C CONTINUE 
C 

c--- USE VALUES FROM KMAT (AT LAST ITERATION A AFTER INVERSICN) TO 
C-- FIND CORRELATIONS WHICH WILL THEN BE SET INTO ARRmY C0RRL<21) IN 
C-- CRDFR ; 2 1 3.1 3,2 4,1 4.2 4,3, ETC. . NOTE THAT FIRST THE 

C-- SOU/O!. ROOTS Or DIAGONAL LLEMENTS WILL OE TAKEN. FOR CONVENIENCE. 
C-- ALSO NOVt THAT ORDER IN THIS CORRELATION ARRAY IS IN TERMS OF 
C-- THE CrOER If. WHICH THE PARAJ'.ETLRS WERE PITIED, NOT THE ORDER IN 
C-- VHIC-I THEY ARE IN A( . ) . 

C 


B6 


DCy9 DO 3001 I»1,NA 

0330 3001 COhr<L(I)«0. 

03 n C 

0362 IO«0 

0363 DO 3005 J-2,NA 

036 » JM=J-1 

0365 DO 3005 I^l.JM 

0336 103lJ«l 

0357 CORRL{10)=XIIAT<I,0)/<XMAT(I,I)*X»1AT(0,0)) 

035S 303D CONTINUE 

0359 C 

0370 C 

0371 C— 3010 REACHED VWEN LINEFIT CONVERGED, PRODUCED PARAMETER ESTIMATES. 

0372 3010 CONTINUE 

0373 C 

0374 C... CHECK LINErIT P,\RAMETERS AGAINST EDIT LIMITS. SIGNAL BV lER.GT.l 
03V5 C 

0376 DO 3012 II-l.NA 

0377 I=OORDR(II) 

0378 VI»A(I) 

0379 3012 IF < {VI.LT,AEDIT<l,l)).0R.(VI.GT.AEDIT(2,l)> ) GO TO 4000 

0380 GO TO 4500 
03S1 C 

0332 4000 lER-l^I 

0333 C 

0334 C--- COMrUTE SWH ESTIMATE, WITH POSSIBILITY OF NEGATIVE SWH . . . 

0J8a 4500 DIFQ=A(3)**2-CALM-*CALM 

0336 Sl7H'rr=0.6*SIGN(SQRT(ABS(DIFQ)),DIFQ) 

0337 C 

0333 RETURN 

0389 END 


B7 




0G31 SUBROUTINE SVMINV (A.N, IFAIL .NROW.P.Q.M) 

0C J2 C. . .PROGRAtV^ED BY: JIM MC1ILLAN - WASC - REVISED 03/0«/78 
£Ci3 C... PURPOSE: TO CWIPUTE THE INVERSE OF A SYTtCTHIC MATRIX 

con DIMENSION A(i\'RO'J,l),P(l),Q(l),M<l) 

GOTO I FAIL = 0 

eOS DO I 1=1, N 

ami 1 M(i) = 1 

mis C... BEG IN CALCULATIO'i. 

0DJ9 DO 14 I = l.N 

C01O C-.-SS.^RCH FOR PIVOT. 


0011 


BIG = 0.0 

C312 


DO 4 J=1,N 

£^13 


TEST = AbS(A(J.J>) 

0014 


IF (TEST-BIG) 4,4,2 

0010 

2 

IF (M(J)> 15,4,3 

0016 

3 

BIG = TEST 

0017 


K « J 

0013 

4 

CONTINUE 

0319 

C... PREPARATION FOR ELIMINATION 

0020 


M(K) = a 

0021 


GiK) = 1.0 / A(K,K) 

0S22 


?<K) = 1.0 

0023 


A(K,K) = 3.0 

L.B2\ 


KPl = K ♦ 1 

JU/'j 


km: = K - 1 

0026 


IF CCMl) 15,8,5 

0027 

5 

DO 7 0=1, KMl 

0023 


0(0) = A(J,K) 

0(32 J 


QCO) = A(J,K) * Q(K) 

0030 


IF (ri(J)) 15,7,6 

OOJl 

6 

G(0) = -0(0) 

00 J 2 

7 

A(0,K) = 0.0 

0033 

8 

IF (K-N) 9,13,15 

0034 

o 

DO 12 0-KPl.N 

0035 


P(0) = A(K,J) 

0036 


I" (N(0)) 15.10,11 

CC37 

13 

P(0) = -?(0) 


11 

0(0) = -A(K,0) • Q(K) 

0039 

12 

A(K,0) = 0.0 


0 *.’,'! C... ELIMINATION FROPEP. 

UL'l! IJ t'O K J=i.:J 

0.541 00 14 K=J,N 

0343 U A(J,K) = A(J,K> ♦ P(J)*Q<IO 

C0J4 RETURN 

0fc4‘> C... ERROR EXIT. 

00tS 15 IFAIL > 1 

e.0'*7 RETURN 

J.'JiJ e:.d 

I 


^ B8 

4 


0 ~J 1 c 

0002 C 

e063 c 
eco4 c 
0CC‘> c 
0C.J5 C 
0037 C 
0338 C 
e.V39 c 

0313 C 
001 1 C 
0312 C 
0013 C 

0314 C 

0015 C 

0016 C 

0017 C 

0018 C 

0019 C 
0320 C 

0021 C 

0022 
0023 
C324 C 
0O2i> 

C32C C 
0027 X 
C32f C 
0025 C 
033.0 
C031 

0032 

0033 
0334 
0035 C 
0.036 
0337 
0:<33 c 
C039 
0340 C 

ca:i c.. 
0 O *>2 c.. 

0043 C.. 

0044 C 
07.;5 c.. 

0045 C.. 

0047 C.. 

0043 C.. 

0049 C.. 

0350 C 

0051 C 

0052 C— 
Ct)'i3 C 
00 54 C 
0055 C 
C3jo X 
C.53V X 
0;j* X 


APPENDIX C. 

SOURCE LISTINGS FOR SUBROUTINES FILLV AND FIELD 

D:^0:GHAVNZ:FILLV3.FR REVISED 10/15/80, APPROX. 1320 HOURS 
G.S.HAVNE WFC/DAS 09/17/80 


THIS VERSION USES THE RADIX 8-4-2 FFT ROUTINES; THE VZRSIOM 
IN FILE FILLV2.fr USES THE RADIX 4-2 RO>JTINES. FILLV3 USES 
THE VALUES OF SVS(514) WHICH ARE PASSED BY A NE'-y LABELLED 
COW10N ARRAY /SVSTTV 

SUBROUTINE EVALUATES SEASAT-1 WAVEFORMS USING FFT TECHNICUES TO 
PERFORM CONVOLUTION OF : 1) SYSTEM POINT-TARGET RESP-ONSE, 

2) SEA-SURFACE ELEVATION DISTRIBUTION, AND 3) FLAT-SEA RESPO^iSE. 
THIS ROUTINE IS SET UP FOR A 512-POIMT TRANSFORM, AND USES FFT 
SUBROUTINES FFA i FFS TOGETHER WITH REQUIRED SUBROUTINES FOR 
FFA AJO FFS (THESE ARE R2TR,R4TR,R8TR,R4SYN,R8SVN, 

ORDl, AND (»D2). 

ALSO USES SUBROUTINES GTSVS, GT5EA2, AND GTFSR2 (A‘4D CFOUT t 
RLOUT ON A /X COMPILE SWITCH). 

SUBROUTINE FILLV(VAL) 

COMPILER STATIC 

INTEGER I 1 ST/0/ . NNP/5 12/.NP2/514/.NC2/257/ 

INTEGER 0DATE(3),JTIME(3) 


REAL VAL(63).SEA(514),FSR(514>,RES(514) 

CCM^iPLEX CTS''S(25' ‘,CTSEA(257),CTFSR(257),CTRES(257),CPHAS 
EQUIVALENCE (S' ‘ ) ,CTSVS( 1 ) ) , (SEAd ) ,CTSEA(1 ) ), (FSR( 1 ) ,CTFSR(1 )) 

1 (RES(1),CTP£S(1)),(A(1),AM?LI),(A(2),TIM0),(A(3),SIGM,4), 

2 (A(4),BGL!N), (A{5),XLMDA),(A(6),XIDEG),(A(7),XKURT) 

CCXX'WN /SSM4N/ A(7),XONST(7),T(63),NA.ITER,SERSQ.CALM,CORRL(2I) 
1 GUESS ( 7 ) .CNSTR ( 7 ) . JCRDR ( 7 ) . AEOIT ( 2,7) 

CCXMON /SrSTry 3VS(514>,NSVS,NSCTR 

. NSYG-^.'' VALUES AT THE POINT AFTER WHICH THE REST OF THE 514 
. SVS(.) VALUES EG’JAL ZERO, AND NSCTR=!N'OEX OF PEAK OF THE 
. POINT-TARGET RtS?C:>iSE. 

W A R N I N G 1 i i 
. SYS(.) CC.'(Ta::,S the TRAI^SFOt'.M after the first call to FILLV3, so 
. 0? NOT change SVS(.) IN THE EXTERNAL PROGRAM. 

. also, if THE INFORMATION IS TO BE PRINTED OLT, DO THIS BEFORE 
. THE FIRST FILLV3 CALL 


NNP IS POINTS IN TRANSF(7Rri (MUST BE POWER OF 2), NP2 IS 2 MORE 
TH.-‘3: NNP, A.'JO NC2 IS » COMPLEX TRAIJSrORM VALUES WHICH WILL RESULT 
FROM THIS. 


CALL DATE(ODATE.I) 

00ATE(3) = CDATE(3) - 1903 
CALL TIKE(OiIf£,I ) 


! 



C939 

X 


VRITE (12,133) 0TINE,O< 


C.“60 

X 

133 

FORMATC FILLV3 DEBUG . 

» 

0£/Sl 

C 




0062 



IF (11ST.NE.0) GO TO 250 


C0u3 



I1S7»1 


0064 



SMSYS=0. 


3035 



DO 150 I«1,NSVS 

t- 

0066 


150 

SMSYS » SMSYS ♦ SYS(I) 


O0S7 



0 ■ NSVS ♦ 1 


>9033 



DO 155 I<>0,514 


0069 

155 

SVS(I) « 0. 


0070 

C 




AT •,I2,2C:',12),* ON * , I2,2< V , 12 ) ) 


• 0371 

0072 

0073 
0374 
0075 
C07C 

0377 

0378 
0079 
0090 
0031 

0332 

0333 
0084 
0335 
038S 
0337 
0388 
0909 
0330 
0-J91 
0392 
0033 
0094 
0395 
0096 
C097 
0038 
0999 
0100 
0101 
0102 

0103 

0104 
Dies 
0106 

0107 

0108 

0109 

0110 
0111 
0112 

0113 

0114 

0115 

0116 
Bin 
r. E 


C IF DESIRED, CAN REPLACE THE BUILT-IN SAMPLED SYSTEM RESPONSE 

C FUNCTION HERE BY ANOTHER ROUTINE GTSYO BY EXECUTING 

C THE FOLLOWING STATEMENT (I.E., REMON t "C“ IN COL.l). 

C CALL STSVSINNP.SYS.SMSYS) 

C 

C— - NOW DO THE TRANSFORM OF SYSTEM IMPULSE RESPONSE 
X WRITE FREE (12) *<15> INPUT REAL SYS(.) :* 

X WRITE FREE (12) * SMSYS » ‘.SMSYS 

X CALL RL0UT(12,SYS,NNP) 

CALL FFA(SVS.NNP) 

WRITE mEE (12) *<15> COMPLEX, TRANSFORMED SYS(.):* 

CALL CPOUT(12,SVS,NC2) 


— N0\^ SET UP SEA SURFACE ELEVATION DISTRIBUTION. 

250 CALL GTSEA(tLN?.SEA.SMSEA) 

C-_- New DO TRANSFORM OF SEA SURFACE ELEVATION DISTRIBUTION. 
X WRITE FREE (12) '<15> INPUT REAL SEA-SURFACE SEA(.): 

X WRITE FREE (12) ’ SMSEA « ',SMSEA 

X CALL RL0UT(12,SEA.NN?) 

CALL FFA(SEA,HNP) 

X WRITE FREE (12) ’<15> COMPLEX, TR.4N3FORMED SEA( . ) : ‘ 

X CALL CPOUT(12,SEA,NC2) 

C 

C 


c--- 

X 

X 

X 


SET UP THE FLAT-SEA RESPONSE FUNCTIOfl. 

CALL GTrSR(NNP,FSR.Si'‘^^S.R) 

TRAi'jSFCRi'I the FLAT-S£A RtSPOiNSE KlWCTION 

WRITE FREE (12) ‘ <m> FLAT-SEA IMPULSE RESPONSE FSR(.):' 
WRITE r.scE (12) ■ SNFSR » ',SMFSR 
CALL RL0,T(12.FSR.r;'JP) 

CALL FFA(FER.f:NR) 

l.'RITE Fi'EE (12) ’<l3> COMPLEX. TRANSFORMED FSR(.):' 

CALL CP0bT(l2,rS.R,NC2) 


500 


-- FORM Amplitude normal izat ion a'wrm, then set up phase wjlti- 
PLIER CE^TA factor D?HI (XNCTR IS FOR POSSIBLE LATER SHIFT OF 
ENTIRE RESULT IF DESIRED. AS NOW SET, FIRST 119 V7J.UES OF RES(.) 
ARE THOSE S=A7.T1INC THE DESIP.ED 63 SAMPLER SET OF SEASAT- 1 AND 
THE CENTER IS AT THE 63TH GATE IF XNOTR IS 0 .) 

XTSCT9 « 0. 

PHI « 0. 

DFHl e -( (XflCTR-24.5-rLOAT(NSCTR) ) ♦ TI.M3/1 . 5625 )• 

1 6.233U53/FLOAT(7tfJP) 

A7KDRM » l./(SM3VS'SKSEA) 

CTHtS(l) « CMPLX(/j\0RM,3. )»CTSVS(:)*CTSEA(t)*CTFSRvl) 

DO 500 le2,NC2 

PHI » PHI * DPHT 

CPHA5 « CMPLX(Ar:CRM’COS(PHl),ANCP:';*SIN(PHI)) 

CTNES(I) = CPHAS«CTSYCll )*CTSEA(I )*CTFSR(I ) 


C2 


0 \- 


V.XLjV' 


0119 

C 

0120 

c 

0121 

C 

0122 

C— 1 

0123 

c 

0124 

X 

0125 

X 

012C> 

0127 

X 

0128 

X 

0129 

X 

0130 

0131 

0132 

c 

0133 

530 

0134 

550 

0135 

C 

0136 

C 

0137 

C 

0138 

0139 
01 if} 
0141 

C 

0142 

610 

0143 

0144 

C 

0145 

620 

0145 

0147 

C 

0143 

C 

0149 

C-— I 

0150 

0151 

0152 

C 


ADOVE WAS A COMPLEX MULTIPLICATION (IN TRANSFCPJ^ DWIAIN) 


WRITE FRCE (12) *<1B> C(3h?LEX, TRANSFORM RESULT RES(.):' 
VttITE FREE (12) ' ANORM ■ ' .ANORM 
CALL CPOOT(12,RES,NC2) 

CALL FES(RcS.NNP) 

V«<ITE FREE (12> '<19> FINAL REAL RESULT RES(.):' 

CALL RLOUT(12,RES,NNP) 

IF (ABS<BSLIN).LE.(AMPLI*0.1E-05)) GO TO B50 
DO 530 I>1,NNP 


TRANSFER DATA TO FINAL OUTPUT ARRAV VAL(63) AND WRITE OUT. 

DO 610 I>1,29 
0»I^I 

VAL(I-»34)-RES(J«61) 

VAL(I)t=RES(J-l) 

DO C20 1-30.34 
VAL(I)bRES(I«28) 


J363 


RETURN 

END 



0001 C DP0:GHAVNE: FI LL03.fr REVISED 09/25/89, APPROX. 093P HOURS 

es'ji c 

CLi.Ta C C.S.HAYNE WFC/DAS 09/25/80 

f^/04 C 

0€3b C F1LLD3 IS EXACTLY SAME AS F1LLD2 EXCEPT FOR THE ADDITION 

C0;76 C OF ANOTHER LA3ELLED COi>lliON AREA /SYSTM/. THIS COMflOr) IS 

Cut37 C FOR COFWUNICATION OF INITIAL SYSTEM POINT-TAP.CET R£Si>CNSE 

00^8 C FUiiCTION VALUES BETWEEN FILLV (FILE FILLV3) AND THE MAIN, 

0009 C CALLING PROGRAM; THIS WAS NECESSARY FOR SSFIT3 AND SIMILAR 

0010 C PROGRAMS... GSH 09/25/80 

0>n 1 C 

0912 C 

0013 C 

0014 C THIS ROUTINE SXME AS FILID.FR EXCEPT FOR DIFFERENT COrG>»N 

0015 C AREA. SSM4N. 

0016 C 

COi: C A SUBROUTINE TO FILL AN ARRAY DRV(63,7) CONTAINING 63 SA1«>LE 

0C1C C VALUES OF THE UP TO 7 DERIVATIVES FOR SEASAT WAVEFORM CASE; 

0019 C THE ORDER OF T.iE DERIVATIVE TERMS IN DRV IS SET BY OORDR(7). 

0020 C FILLD IS A COMPANION TO FILLV; IT REQUIRES iHAT FILLV HAS BEEN 

0021 C CALLED ALREADY. SEE FURTHER COMfiaJTS IN THE FITTING ROUTINE 

0022 C SWhFIT WHICH CALLS FILLD AND FILLV. THESE ROUTINES ARE PART 

0023 C OF THE NEW PROGRAMS USING THE FFT-CON^VOLUTION TECKYIQUE WITHIN 

0024 C WAVEFORM GENERATiai. DEBUG PRINTOL'T TO (ASSUMED ALREADY OPENED) 

0925 C DEVICE 12 IS PROVIDED aN X COMPILE SWITCH. 

0926 C 
0327 C 

0S2C C FILLD2 REQUIRES SUBROUTINES FILLV, GTSYS,CTSEA2,GTFSR2,FR27R,FR4T 

0029 C FR4SVN,F0.TD1 , AND f0R02. ' 

f;>30 C 

£031 C THE DERIVATIVES ARE NUMERICALLY FOUND BY STEPPING A(7), THE 

a~Z2 C PAT. f METER VALUES BY AN AMOUNT SET BY STPRM (EXCEPT THAT FIRST 

0033 C AND FOURTH TERMS. AMPLITUDE 4 BASELINE, ARE FOUND DIFFERENTLY). 

C"3 ) C 

0UJ5 SUDRCUTIMi FiLLD(VAL .DRV. STPRM) 

0936 C0:PILER STATIC 

r:37 c 

00i3 COMMON /SYSTM/ SVS(514) .NSYS.NSCTR 

0939 C 

0940 CDMMOiV /SSM4N/ A( 7 » ,Xa!ST(7 ) ,T(63 ) ,tO\, ITER.SERSQ.CALM.CORRLIEl ) , 

0941 1 CUESS(7).C'.STR(7),UCRD?.(7).AEOIT(2,7) 

09 42 C 

UC‘*i EC'JIVALEMCE (A.*-:^LI , A( 1 ) ) , (DSLIN, A(4) ) 

0944 C 
0745 C 

0346 REAL VAL(C3> .CRV(63,7) ,STPRM(7) 

03 »7 X PEAL XT.'’J>(fe3> 

0948 C 
B'JiO C 

CO60 X WRITE FREE (12) • DEBUG IN FILLD; CALL FILLV A PRINT RESULTS 

0051 X CALL FILLVCU.O.) 

0952 X CALL RLO'.iT(12.VAL.63) 

0.;S3 C 
0954 C 

9Dab C— LOO? TO 40C FOR THE NA DERIVATIVES NtCDEO. 

0*156 C 

0357 CO 400 K-l.NA 

C/E? OrOORDR(K) 


C4 


C 

0363 C- — 
0061 C- — 
03-1 2 C 

0063 X 

0064 C 
0365 

036 S C 

0067 C— 
3363 100 

0363 X 
3370 

0371 

0ii72 

0373 110 

0374 
037' C 
0376 C— 
03/7 200 

0078 X 

0079 

0633 210 

0331 
03S2 C 
0<JO3 C™ 
0034 300 

0033 X 
03J6 
0337 
0033 

001' 3 C 

0093 C-” 
0391 X 
0032 

0393 310 

0094 C 

0395 C--- 

Ovse c 

0037 X 
0090 C 
0059 
DliO X 
0101 X 
01f/> C 
0103 C 
010'. C— 

0105 C 

0106 
010V 
0108 

0109 X 

0110 323 

0111 C 

0112 
0in C 

0114 C 

0115 X 

one X 

0117 C 

OllP c 


DO NUf.L*RICAL STEP PROCEOURc AT 320 FOR EVERYTHING EXCEPT AMPLI AND 
BSLIN; FO% THESE UE KNOW EASIER WAV TO CET ANSWER. 

WRITE FREE (12) * FIELD DEBUG: AT COMPUTED CO TO, 0 * K « ‘.J.IC 

GO TO (103. 300. 300, 200, 300, 3J0. 300), 0 

100 IS STATErXNT REACHED WHEN AMPLI OEP.IV IS NEEDED. 

CONTINUE 

WRITE FREE (12) • FIELD DEBUG: AT STATEMENT * 100' 

TMP^AMPLI 

IF (TMP.LT.I.E-05) TMP=l.E-05 
DO 110 I>1,63 

DRV ( 1 ,K ) ■ ( VAL ( I ) -BSL IN )/TMP 
GO TO 400 

200 IS STATEMENT REACHED 1/HEN BSLIN OERIV IS NEEDED. 

CONTINUE 

WRITE FREE(12) ‘ FIELD DE&'JG: AT STATEMENT » 20C' 

DO 210 i>l,63 
DRV(I,K)«1.0 
GO TO 430 

303 IS STATEMENT FOR NUMERICAL STEP ESTTMATION OF THE DERIVATIVE. 
CONTINUE 

WRITE FREE (12) * FIELD OEB'JG: AT STATEMENT • 3C3‘ 

ATMR=A(0) 

STEP»STPRM(0) 

A(J)»A(0)*STE? 

ABOVE STATEI»£NTS SWAPPED IN NEW VALUE FOR THE DESIRED A(.) 

NOW STORE ELEMENTS OF VAL(.) TEMPORARILY IN DRV(.,J) 

WRITE FREE (12) * FIELD DEGUG: 0 « A(J) ■ '.J,A(J) 

DO 310 I>1,63 
DRV(1.K)»VAL(I) 

CET NEW WAVEFORM WITH MEW (TEM‘ORARV, STEPPED) PA.RAMTTER VALUES. 
'■/RITE FREE (12) ' <LF> TEMPORARY A( . ) :‘,(A(I), I»l,7) 

CALL FlLLV(VAi) 

WRITE Ff.EE (12) * AND FOR THESE VALUES, FUNCTION IS :' 

CALL RLOOT(:2.VAL.S3) 


NOW FCRM DEPIV. ESTIMATES & S'VAP VALUES BACK TO ORIGINAL.. 

DO 330 I«1.63 
Tr'rf>=CR'/(I ,K) 

DI:V(1,K)-(VAL(1)-TM?)/STEP 

XTMP(1)«DRV(I,K) 

VAL(I >-'TMP 


A( J)«ATf3» 


WRITE FREE (12) ’ <LF> STEP-DERIV. VALUES ARE 
CALL RLOJT(12,XTMP,63) 


ffll9 

X 

WRITE FREE (12) ' REST0.5ED FUNCTION VALUES 

0120 

X 

CALL RLOUT(12,VAL,63) 

0121 

C 


01 2r 

C 


0123 

C— - 

BOTTOM 0" f’AIN DERIVATIVE LOOP 

0124 

c 


012S 

400 

CONTINUE 

0126 

c 


0127 


RETURN 

0128 


END 


APPENDIX D 


SOURCE LISTINGS FOR SUBROUTINES GTFSk AND GTSEA 


0S31 

0:752 

00'yj 

eoDX 

0053 

0056 

6007 

0008 

00.79 

0010 

0011 

0012 

0013 

0014 
6015 
0O1C 
0617 
6018 

001 5 
0020 
0521 
(.072 
G'720 
0024 
6375 
6026 
P077 
6023 
0-779 
O'ljO 
0031 
6032 
0033 
0534 

0035 

0036 
0J37 
0033 
0CJ9 
0010 
6541 
6042 
0J13 
0044 
OOC5 
00 ;6 
6047 
0643 
0049 
00.40 
6031 
C0J2 
6653 
6l>34 
o;55 
6036 
6 .'a7 
6' :3 


C DP6:GHAYNL: GTFSR2.fr REVISED 11/13/80, 1505 HRS 

C 

C C.S.HAVKE APPL. SCI. ASSOC *S 07/16/85 

C 

C THIS VERSION ADMITS A NEGATIVE ANCLE, BUT FAKES THE RESULT 

C TO GIVE A FASTER -THAN-ZERO-DECREES PLATEAU DECAY. 07/U/80 

C 

C SUBROUTINE 6TFSR2 FILLS THE FLAT-SURFACE RESPONSE, USING 10 

C TERM ONLY FROM EXPANSION IN GARY B.ROUN'S PAPER. A POWER SERIES 

C FROM ABRAMOWITZ « STEGUN IS USED TO EVALUATE 10. THIS VERSION OF 

C GTFSR USES 230 NON-ZERO VALUES OF THE FLAT-SURFACE RESPONSE, AMO 
C ASSUMES THAT NNP > 230. GTFSR2 IS DERIVED FROM GTFSR BUT IS 

C FOR THE NEW FFT-BASED CONVOLUTION USING ONLY 5I2-PCINT 

C FFT. 

C 

SUBROUTINE GTFSR ( NNP, FSR.SMrSP) 

REAL FSR(l) 

C0«40f4 /SSK4N/ AMPLI ,TlM0.SIG:‘1A,3SLIN,XLrt)A,XIDEC. VKURT 
NFSR - 230 
SMF3R « 0. 

OT « 1.5625 
T « -DT/2. 

C 

C... TEST FOR (IMPOSSIBLE) NEGATIVE ANGLE; IF PRESENT, CHOOSE BRANCH WHICH 
C... EFFECTIVELV INCREASES THE OLTA AT < ZERO DEGREES POINTING (CHANGES 
C. . . M4DE ON 07/16/80) 

C 

IF (XIDEC.CT.O. ) CO TO 20 
C 

C... THE FOLLOWING STATEMENT CAUSES DLTA TO INCREASE BV A FACTOR OF 
C.. TWO FCn 1 DEGREE (FICTICIOUS) NEGATIVE SEASAT ANCLE... 

DLTA ■ 2.66496E-03"(I.-X1DEG) 

DO 15 Jal.NFSR 

T ■ T ♦ OT 

Z » AM?LI»EXP(-DLTA‘T) 

FSR(O) » Z 

15 SMFSR -s SMrSR ♦ Z 
GO TO 35 
C 

20 X2RAD • X:0E3/28.6«789 

BtTA ■ 4.35331«S:N(X2RAD) 

DLTA ■ 2.SS4.6E-03*CO3(X2RAD) 

DO 30 3«1,\-SR 


T ■ T • DT 
Z • OLTA-SCRT(T) 

C SELECT WHICH Or THE TWO SERIES TO USE FCR 10(Z) 

IF (Z.GT.3.75) CO TO 23 
Z » Z-r/14.0626 

C 3.75«*2 ■ 14.6625 

A ■ 1. ♦ Z«(3.515623*Z»(3.039942*Z“(1. 206749 
1 *Z*(0.?659732*Z*(0.0360763*Z*0.064S813>)))) 

CO TO 27 

23 A ■ F.X‘*(Z)/SCRT(Z) 

2 ■ J.75/Z 

A ■ A*(0.3989423-Z*(0.03983024«Z*(0. 00362018 

1 -Z* (0. f)?l 63801 -Z*(0. 0103 1555-Z*(0.02282967-Z* (0.028958 12 

2 -Z‘ (0.01 787654-2«0, 03426059) ))))))) 

77 Z ■ AliPLI*EXP(-DLTA'‘T)‘A 


lK50 

0aiii 

ca.2 

cr.(i3 

0CS4 

0iu& 

C0C6 

0Oi7 

CJ63 

eei9 

0070 

6371 

0672 

0373 


FSR(J) ■ 2 
30 STFSH ■ ShiFS.'i ♦ 2 
C 
C 

35 CONTINUE 
C 

C FILL REST OF THE AK71XV WITH ZEROES 

K ■ NFSR ♦ 1 
NP2 ■ NNP ♦ 2 
DO 40 0 ■ K.NP2 
49 FSR(J) ■ 0. 

IlST « 1 
RETURN 
C 

END 


02 


REVISED 10/15/10, 1400 HRS 


0C01 

0Cn?. 

Ci.OZ 
LC94 
0039 
0C7i 
BjZ7 
0J08 
00 V9 
C310 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 
0019 
0320 
0021 
0322 

0023 

0024 
0026 
eC2i 
01*27 
0028 
CL29 
0030 
P031 
0032 
DC33 
00 J 4 
0035 
CJ16 
0337 
0038 
0339 
Car.3 

0041 

0042 
O.J43 
e.*44 
CO 1b 
0016 
0017 

eo-.« 

C J >9 
C’-m 
0531 
rj52 
0.''33 
0:54 
Ciab 
OJi.6 
C.'jI 

o-'j'; 


C DPO:CHAVNE: CTSCA2.fr 

C 

C G.S.HAVNC APPL. SCI. ASSOC 'S 06/20/80 

C 

C DISCOVERED IMPORTANT SIGN ERROR IN THIS ROUTINE ON 10/15/80; 

C CH.\NCED IT AT ABOUT 1400 HOURS. PRIOR TO THIS. ALL 

C CTSEA2 RESULTS WERE RETURNING OPPOSITE SIGN OF SKEWIVESS TO 

C THAT INTEfOEO... 

C 

C SUBROUTINE GTSEA2 FILLS ARRAY WITH A SKEWED GAUSSIAN SURFACE 

C ELEVATION DISTRIBUTION. CENTERED ON THE SAMPLE NUTOER 86. 

C ZEROES ENTERED IN ALL OTHER ELEMENTS THAN IN INTERVAL 1 - 171. 

C ASSUMES NNP > 171 . THIS ROUTINE SAME AS CTSEA.FR EXCEPT FOR DIFFERENT 

C COt-'.'lON NAME AND FOR THE DISTRIBUTION CENTERED AT 171. IT IS INTT/UPEP 

C FOR USE IN THE 512-POINT FFT PROCESSES. 

C 

SUBROUTINE GTScAvNNP.SEA.SMSEA) 
real SEA(2) 

COMMON /SSM4N/ AMPLI ,TIM3,SIG7U.BSLIN.XLFLA.XIDEG.XKURT 

C- CONVERT SEA SIGMA TO » GATE INTERVALS 

XNGTS - SIGMA/1.5625 

C zero width IS NOT ALLOWED 

IF (XNGTS.lt. 0.001) XNGTS ■ 0.001 
WGTS ■ 0. 

C ESTABLISH CENTER AT 86TH CATE 

NCTR ■ 86 
SEA(fiCTR) ■ 1. 

SMSEA « 1. 


K ■ NCTR - 1 
X5 ■ XLFOA/6. 

XflDX ■ 0. 

C FILL NCN-ZEPO ELEMENTS OF THE ARRAY 

DO 20 J>1.K 

XrJOX > XNDX * 1. 

WGTS « XNDX/XNGTS 
2 ■ -WCTS*WGTS/2. 

IF (2.LT.-8C. ) CO TO 10 
A ■ EXP(Z) 

2 ■ X6»WSTS“(WCTS*UCTS-3.) 

C 

C 

C... REVERSED SIG7S IN A1.A2 BELOW AFTER FINDING SIGN ERROR ON 
C... 10/15/80 

C 

A1 ■ A*<1. - Z) 

A2 ■ A-(l. • 2) 

C 

C 

C 



CO 

TO 15 



10 

A1 

■ 0. 




A2 

•> 0. 



15 

SEA(NCTR- 

J) 

« A1 


SEACNCTR* 

0) 

■ A2 

20 

SMSLA 

• SMSEA * A1 


C NOW FILL ZERO ELF71ENTS 0- THE ARRAY 

r » 2-NClR 
NP2 • NN? ♦ 2 


C'3? 00 30 0-K.NP2 

OCiC 30 SLA(J) > 0. 
r>961 KETURU 

COL 2 C 


APPENDIX E 


SOURCE LISTINGS FOR SUBROUTINES FFA AND FFS 


eaOl C OP0:CHAVNE:FFA.FP REVISCU 07/03/S0, APPROX. 10SB HOURS 

ei’Ot c 

C G.S.HAVNE APPL .SCI .AS >OC‘ S. 0//01/S0 

C 

C'jab C THE FOLLOWING IS ONE OF THE SET OF ROUTINES FOR FAST FOURIER 

rOCb C TR/MSFORM OF REAL DATA SEQUENCE AS DESCRIOEO IN SECTION 1.2 OF 

enn c "programs for oighal sic.nal processing." ed. bv digital signal 

0008 C PROCESSING CCMMITTEE OF THE IEEE ASSP, PUDLISHEO BY IEEE PRESS, 

0009 C NY, 1979. THESE ROUTINES ARE COLLECTIVELY THE "FFA-FFS PACKACE" 

W\0 C WHICH INCLUDES: FFA. FFS. R2TR, R4TR, R8TR, R4SYN, R8SYN, ORDl , 

0:91 1 C AND 0KD2. 

0012 C 

0013 C - 

03U C SUBROUTINE: FFA 

0019 C FAST FOURIER ANALYSIS SUBROUTINE 

001S c 

CD17 C 

0018 SUBROUTINE FFA(B, NFFT) 

0019 C 

C02O C THIS SUaSOL'TINE REPLACES THE REAL VECTOR B(K>, (K»l,2 N). 

0021 C WITH ITS FINITE DISCRETE FOURIER TRANSFORM. THE DC TERM IS 

0022 C RETUR/.'L'D IN LOCATiai B(l) WITH B<2) SET TO 0. THEREAFTER, THE 
0323 C JTH HAT.KONIC IS REfURNED AS A COMPLEX NUMBER STORED AS 

0024 C B(2*0*l) ♦ I B(2*0*2). NOTE THAT THE N/2 H/J?MONir IS RETURNED 

0025 C IN B(N*1) WITH 3(M*2) SET TO 0. HENCE. B MUST BE DIMENSIONED 

0326 C TO SI2E N»2. 

0027 C SUBROUIINE IS CALLED AS FFA (B.N) WHERE N>2**M AND B IS AN 

0020 C K TEk.i P.t.«L ARRAY. A REAL-VALUED, R/>DIX 8 AiGCRITH.H IS USED 

0029 C WITH IN-PLACE REORDERING A\’D THE TRIG FU7.CTI0NS ARE CCM^UTcO AS 
£330 C NEEDED 

C031 C 

0032 DIMENSION B(2) 

0033 CO«,'»N /CON/ PH. ?7 , P7T'.X), 022. S22. P12 
033.. C 

0035 C IW IS A MACHINE-OEPEr.'DENT URITE-DEVICE NU.rSCR 

C936 C 

0037 1V;»12 

0038 C 

C.T39 P?I ■ 4."AT.jm. ) 

C3/-0 Pie » PI 1/9. 

con P 7 « i./scs:t(2.) 

DC'.Z P7TW0 ■ 2."»7 

0743 C22 ■ C3S:?;e) 

c’o:; S22 ■ c!N<?:8) 

OC45 P12 » :."?II 

0JM6 N • : 

F347 DO 10 I«1.10 

0043 M « I 

0.7*9 N ■ N*2 

0:.5D IF (N.EO.NFrT) CO TO 20 

0051 10 CONTIf.ME 

Ct352 W?17E (IW.9999) 

C063 9999 FOR^HaTC NFFT NOT A POWER OF 2 FOR FFA*) 

0Jj4 STOP 

ecus 20 CUNTINUC 

035b N8P03 » M/3 

03a7 C 

.1 C DO A RaDX 2 OR RADIX 4 ITERATION; FIRST IF C:i* IS FEOUIRED 


cr>) 

C 



0PiO 



IF <M-NCPOW*3-l) 59, 40, 30 


30 


NN «: 4 

CCi-: 



I NT 3 N/NN 




CALL R4TR(INT, BU), B(1NT*1), P(2»IiMT*l), B(3*IN*>1)) 

fj.-j '.1 4 



GO TO 60 

£l/bO 

40 


r.TJ = 2 

0CG3 



INT a N/NN 

WO 7 



CALL R2TR(INT, B(l), 3(INT*D) 

0uS2 



Uu TO 60 

POC9 

50 


NN a 1 

0J7O 

C 



UJ71 

C PERFORM P^IX 8 ITERATIONS 

0072 

C 



0373 

60 


IF (N8POW) 90, 90. 70 

0374 

70 


DO 80 ITal.N8POW 

0370 



NN a NN«8 

0376 



INT a N/NN 

03/7 



CALL R8TR(INT,NH.B(1),B(INT»1),B<2*1NT^1),B(3*INT*1) 

0078 


1 

B(4-:NT*1),B(5*IN1>1),B(6»INT*1),B<7-INT*1),B(1>, 

0379 


2 

B(INT*1) ,B(2*INT*1),B(3«INT*1),B(4*1NT»1),B(5*INT 

0330 


3 

B(6‘INT*1),B{7*INT*1)) 

0981 

30 


CONTINUE 

0332 

C 



0333 

C PERFORM IN-PLACE REORDERING 

0334 

C 



0333 

90 


CALL ORDKM, B) 

0036 



CALL ORD2(M, B) 

C587 



T a B(2) 

0088 



B(2) a 0. 

0039 



B(NFFT*1) a T 

03S0 



b<NFFT^2) a 0, 

0391 



DO ice I=4.NFrT,2 

0332 



B(I) a -3(1) 

9093 

103 


CONTINUE 

0394 



RcfURN 

0090 



END 


tJRIGINAL PAGE IS 
OF P(X)R QUALITY 


0001 

0217. 

0-333 

0034 

iKJ& 

0006 

0607 

0SO8 

0039 

0010 

001! 

U012 

6913 

0014 

0015 

0016 
0017 
0313 
0119 
0020 
0021 
0022 

0023 

0024 
0325 
0026 
0127 
0023 
0029 

0130 

0131 

0032 

0033 
OU34 
0 335 
C 336 

033 7 
0-333 
003 7 

0-0 to 

COM 

0J»2 

00'? 

0044 

00 46 
Ot W 

aon 

0049 

0050 
0-031 
2052 
01*13 

0034 
0-:35 
0556 
OJj? 

1 . ViJ 


C Dr*0:GHAVNE:FKS.FR REVISED 37/03/SO, APPROX. 1130 HOURS 

C 

C G-S-HAWIE APPL.SCI.ASSOC'S. 07/02/90 

C 

C THE FOLLOWING IS ONE OF THE SET OF ROUTINES FCR FAST FOURIER 

C TRANSFORM OF REAL DATA SEQUENCl AS DESCI^IDED IN SECTICN 1.2 OF 

C “;»ROGRAMS FOR DIGITAL SIGNAL PROCESSING," ED. BY DIGITAL SiGNAL 

C PROCESSING COrniTTEE OK THE IEEE ASSP, PUBLISHED BY IEEE Pl-^tSS, 

C NY. 1979. THESE ROUTINES COvLECTIVELV ARE "KFA-FFS PACKAGE* 

C WHICH INCLUDES: FFA, FFS. R2TR, R4TR, R8TR, R4SVN, HUFVN, OKDl, 

C /V4D 0RD2. 

C 

C 

C 

C SUBROUTINE: FFS 

C FAST FOURIER SYNTHESIS SUBROUTINE 
C R/OIX 8-4-2 

C 

C 

SUBROUTINE FFSlB, NFFT) 

C 

C THIS SUBROUTINE SYNTHESIZES THE REAL VECTOR B(K), WHERE 

C K-1.2 N. THE INITIAL KOLHIER COEFFICIENTS ARE PLACED IN 

C THE B ARRAY OF SIZE N*2. THE DC TERM IS IN B(l) WITH 
C D(2) EQUAL TO 0. 

C THE OTH HARMCNIC IS STORED AS B(2»0*l) ♦ I B(2*0*2). 

C THE N/2 HARMONIC IS IN B(M*1) WITH |\(N*X) EQUAL TO 0. 

C THE SUOROUTINE IS CALLED AS FFS(B.N) \WERE N'2**M AND 
C D iS THE N TERM REAL ARRAY DISCUSSED /30VC. 

C 


IW 


DIMENSION P(2) 

CO.mN /CONI/ PIl, P7. P7T»Y0, C22. S22, PK 
IS A MACHIiNE-DEPET.’DcNT i.’RnE-DcVK,'' NO^SER 
12 • 12 


10 

1599 

20 

30 


PII ■ 4.*ATAN(1.) 

PIP » PlI/3. 

P; « 1 ./SGST<2. ) 

P7FW0 » 2.*P7 
C22 ‘ COa(PIS) 

S22 ■ SIN p:j> 

?!2 » 2.*?i: 

N ■ 1 

DO le i«:.!5 
M • ! 

N » N'2 

if (N.EQ.N’FT) GO TO 20 

CONTINUE 
WRITE (1W.9999> 

fCR.lUC NFFT NOT A POWER O 2 rOR FFS') 
STOP 

CC.NTlNUe 

B(?) ■ P(NFFT-l) 

DO 3C !■! .NFFT 

Btl ) « D(I )/rLau(Nfrn 

ccnvic'JE 






C3 


0». oO 

GC51 

cc::2 

J 

f.Oi ’ 
r '.'j 
UyiH 
CJsi7 
■JOjS 
OC^SO 
D07:J 
i7071 
OG72 

0073 

0074 
00 70 
0076 
60 V 7 
0373 
0379 
03J0 
0031 
003’ 
Ov/33 
00S4 
C*.30 
Oh/3G 
6Cit7 
0088 
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DO 4<5 I = 4.NKFT,2 

a<l) » -B(l) 

40 CONiriNUE 

K'S»>OW « .1/3 
C 

C RECROCR THE INPUT FOURIER COtlFF ICIENTS 
C 


C 

C 

C 

C 


CALL 0RD2(M, B) 

CALL ORDKM, B) 

IF (N8POW.EQ.0) GO TO 00 

PERFORM THE RADIX 0 ITERATIONS 

MM » N 

DO 50 IT-1, N3P0'V 

I NT « N/NN 

CALL R8SVN<1NT,NN,B,D{1NT^1),B(2*INT*1),B(3»IMT^I), 

1 B.4*INT*1),B(5»INT*1),B(6*INT*1),B(7«INT^1),B<1), 

2 EaNT*l).B(2*INT*l),B<3*INT^l),B(4«lMr*l).B(5*lNT^l), 

3 B(6»INT*1).B(7*INT*1)) 

NN = NN/8 


50 CONTINUE 
C 

C DO A RADIX 2 OR RADIX 4 ITERATION IF ONE IS REQUIRED 
C 


60 IF (M-N8POW-3-1 > 90, 80, 70 

70 INT * N/4 

CALL R4SVNC1NT, C(l), BUNT*!), B<2*INT*1), B(3*INT*l>) 
GO TO 90 
80 INT - N/2 

CALL R2TR(INT. B(l), B(1NT*1)) 

90 RETURN 

END 
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